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Abstract 

Current  flight  reference  systems  are  vulnerable  to  GPS  jamming  and  also  lack 
the  accuracy  required  to  test  new  systems.  Psendolites  can  augment  flight  reference 
systems  by  improving  accuracy,  especially  in  the  presence  of  GPS  jamming.  This 
thesis  evaluates  a  pseudolite-based  flight  reference  system  which  applies  and  adapts 
carrier-phase  differential  GPS  techniques.  The  algorithm  developed  in  this  thesis 
utilizes  an  extended  Kalman  filter  along  with  carrier-phase  ambiguity  resolution 
techniques. 

A  simulation  of  the  pseudolite-based  positioning  system  realistically  models 
measurement  noise,  multipath,  pseudolite  position  errors,  and  tropospheric  delay. 
A  comparative  evaluation  of  the  algorithms  performance  for  single  and  widelane 
frequency  measurements  is  conducted  in  addition  to  a  sensitivity  analysis  for  each 
measurement  error  source,  in  order  to  determine  design  tradeoffs.  Other  analyses 
included  the  use  of  optimal  smoothing,  non-linear  filtering  techniques,  and  code 
averaging.  Specific  emphasis  is  given  to  two  alternate  methods,  both  developed  in 
this  research,  for  handling  the  residual  tropospheric  error  after  applying  a  standard 
tropospheric  model. 

Results  indicate  that  the  algorithm  is  capable  of  accurately  resolving  the  pseu¬ 
dolite  carrier-phase  ambiguities,  and  providing  a  highly  accurate  (centimeter-level) 
navigation  solution.  The  filter  enhancements,  particularly  the  optimal  smoother  and 
tropospheric  error  reduction  methods,  improved  filter  performance  significantly. 


DEVELOPMENT  AND  SIMULATION  OF  A 


PSEUDOLITE-BASED  FLIGHT  REFERENCE  SYSTEM 


I.  Introduction 


1.1  Background, 

Applications  for  the  Global  Positioning  System  (GPS)  have  increased  tremen¬ 
dously  since  its  inception,  including  the  development  of  many  differential  GPS  (DGPS) 
techniques.  Differential  GPS  performs  relative  positioning  between  two  or  more  re¬ 
ceivers  by  calculating  and  removing  sources  of  errors  that  are  common  between  re¬ 
ceivers.  The  integration  of  a  GPS  receiver  and  an  Inertial  Navigation  System  (INS) 
is  another  application  that  has  produced  accurate  and  robust  navigation  systems. 

One  of  the  most  advanced  navigation  systems  is  the  modern  flight  reference 
system  operated  by  the  746th  Test  Squadron,  Holloman  AFB,  NM,  which  is  used 
to  test  and  evaluate  new  flight  navigation  systems.  To  be  useful,  a  flight  reference 
system  should  have  an  order  of  magnitude  greater  accuracy  than  the  system  under 
test,  because  the  output  from  the  reference  system  is  regarded  as  the  truth.  Any 
degradation  in  reference  system  performance  will  invalidate  the  evaluation  of  the 
system  under  test.  The  flight  reference  system  has  evolved  through  the  years  from 
radar  tracking,  ground-based  camera  and  aircraft  transponders,  to  the  current  sys¬ 
tem  of  DGPS  integrated  with  an  inertial  unit,  barometric  altimeter,  and  a  ground 
transponder /interrogator  system  [14].  The  current  reference  system  used  by  the 
746th  Test  Squadron’s  Central  Inertial  Guidance  Test  Facility  (CIGTF)  is  called  the 
CIGTF  High  Accuracy  Post-processing  Reference  System  (CHAPS)  [17]. 
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The  type  of  differential  GPS  that  CHAPS  uses  is  called  carrier-phase  DGPS, 
which  relies  on  measuring  the  carrier  component  of  GPS  signals.  The  carrier-phase 
signal  provides  a  relative  measurement  because  the  total  number  of  carrier  phase 
cycles  are  not  known.  Phase  ambiguities  exist,  which  are  the  unknown  number  of 
carrier-cycles  present  at  the  start  of  the  signal  integration  [28].  The  carrier  signal 
can  be  broken  up  into  two  segments  that  are  separated  by  the  point  in  time  that 
signal  integration  starts.  The  first  segment  consists  of  the  unknown  integer  number 
of  cycles  up  to  the  point  of  signal  integration.  The  second  segment  consists  of  the 
measured  carrier-phase,  which  is  not  constrained  to  be  an  integer.  The  highest  level 
of  accuracy  is  attained  when  the  unknown  number  of  cycles  before  signal  integration 
is  determined. 

Carrier-phase  DGPS  can  be  categorized  into  two  classes  based  on  how  the 
estimation  of  the  unknown  integer  cycles  is  performed.  Floating-point  carrier-phase 
techniques  estimate  the  integer  number  of  cycles  as  floating-point  numbers,  without 
forming  integer  ambiguity  values.  Fixed-integer  carrier-phase  techniques  select  a 
set  of  integer  ambiguities  through  a  process  called  ambiguity  resolution  [15].  If  the 
correct  ambiguities  are  selected,  a  fixed-integer  solution  achieves  greater  accuracy 
than  floating-point  solution.  Fixed-integer  solutions  are  vulnerable  to  selecting  the 
incorrect  integers,  which  result  in  degraded  performance. 

CHAPS  currently  faces  two  challenges:  accuracy  during  periods  of  GPS  jam¬ 
ming  and  accuracy  when  a  GPS  signal  is  available.  Operation  in  the  presence  of 
GPS  signal  interference  impedes  CHAPS  from  using  its  most  accurate  sensor.  When 
jamming  denies  CHAPS  from  using  GPS  measurements,  it  must  rely  on  its  other  sen¬ 
sors,  primarily  the  INS.  INS  accuracy  degrades  over  time,  causing  the  performance 
of  CHAPS  to  suffer.  Although  post-processing  techniques  are  applied  to  reduce  the 
impact  of  INS  errors,  CHAPS  cannot  maintain  centimeter  level  accuracy  during  pe¬ 
riods  of  GPS  jamming.  The  second  challenge  facing  CHAPS  is  accuracy  when  GPS 
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is  available.  As  new  systems  become  more  accurate,  CHAPS  must  also  improve  its 
accuracy  to  be  a  useful  reference. 

One  technology  that  potentially  solves  the  challenges  of  reference  system  accu¬ 
racy  is  the  use  of  pseudolites  (or  pseudo-satellites)  [25] .  Pseudolites  are  ground-based 
transmitters  that  send  GPS-like  signals  which  can  be  received  with  GPS  receivers 
that  are  adapted  for  pseudolite  signals.  DGPS  techniques,  such  as  carrier-phase  am¬ 
biguity  resolution,  can  be  adapted  and  applied  for  pseudolite  navigation.  Pseudolites 
have  the  flexibility  of  operating  at  various  frequencies,  which  allow  them  to  avoid 
GPS  jamming  signals. 

The  aiding  with  pseudolites  will  increase  the  accuracy  of  CHAPS  or  other 
flight  reference  systems  when  GPS  jamming  signals  are  present,  and  also  during 
periods  of  normal  GPS  operation.  Pseudolites  potentially  can  provide  CHAPS  with 
a  navigation  source  that  maintains  centimeter-level  positioning  accuracies  during 
periods  of  GPS  jamming.  The  increase  in  accuracy  during  normal  (non-jamming) 
GPS  operation  is  the  result  of  CHAPS  having  access  to  two  highly  accurate  sensors, 
as  compared  to  just  one  when  pseudolites  are  not  used.  A  system  that  uses  two 
sensors  with  roughly  the  same  accuracy  can  expect  to  see  a  1/ y/2  factor  improvement 
in  accuracy  over  just  using  one  sensor.  That  is  nearly  a  30  percent  improvement, 
assuming  that  both  sensors  are  independent.  The  errors  between  GPS  and  pseudolite 
signals  are  not  completely  independent,  but  a  practical  system  would  still  show 
improvement  over  DGPS-only  navigation. 

1.2  Problem.  Definition 

The  primary  goal  of  this  thesis  is  the  development  and  testing  of  an  algo¬ 
rithm  that  resolves  the  carrier-phase  ambiguities  of  ground-based  GPS  transmitters 
called  pseudolites.  Figure  1.1  depicts  the  process  of  resolving  the  ambiguities  first  to 
floating-point  numbers  in  a  Kalman  filter,  resolving  the  phase  ambiguities  to  integer 
numbers,  and  then  passing  the  carrier-phase  measurements  with  resolved  ambigui- 
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ties  to  aid  the  flight  reference  system.  This  research  includes  the  simulation  of  the 
pseudolite  environment,  along  with  the  creation  of  realistic  errors  in  the  pseudolite 
code  and  phase  measurements.  The  secondary  goals  include  analysis  of  the  effect  of 
each  error  source,  and  the  application  of  new  methodologies  for  dealing  with  pseu¬ 
dolite  error  sources.  The  objective  of  this  work  is  to  develop  an  algorithm  that  can 
improve  both  the  accuracy  and  robustness  to  jamming  of  a  flight  reference  system 


Floating  point 
ambiguities  and 


Pseudolite  carrier- 
phase  with  resolved 
ambiguities 


Flight  “truth” 


Figure  1.1  Pseudolite  Aided  Flight  Reference  System 

1.3  Related  Research 

Raquet  et  al.  [25]  conducted  an  early  test  of  a  pseudolite-only  flight  reference 
system.  This  work  was  accomplished  at  Holloman  AFB  under  the  partnership  of  the 
746th  Test  Squadron,  NovAtel  Communications,  Stanford  Telecom,  and  the  Univer¬ 
sity  of  Calgary.  This  proof  of  concept  involved  an  ’’inverted”  mode  of  operation  in 
which  the  position  of  the  pseudolite  is  solved  in  relation  to  an  array  of  stationary 
receivers  that  are  placed  at  known  locations. 
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NovAtel  Communications  and  Stanford  Telecom  continued  work  with  pseudo- 
lite  navigation  by  conducting  a  follow-on  test  to  duplicate  some  of  the  results  from 
the  Holloman  proof-of-concept  test.  They  prototyped  a  GPS/pseudolite  system  that 
allowed  coverage  during  times  of  reduced  GPS  availability  [12]. 

Although  pseudolite  signals  are  very  similar  to  GPS  signals,  many  assumptions 
that  are  made  for  GPS  navigation  cannot  be  applied  to  pseudolite  operations.  Sec¬ 
tion  2.6.1  details  the  differences  between  GPS  and  pseudolite  systems.  Dai  et  al.  [8] 
addressed  some  of  the  challenges  that  pseudolites  present  by  developing  unique  mod¬ 
elling  strategies  to  deal  with  pseudolite  error  sources.  They  also  analyzed  the  impact 
of  pseudolite-user  geometry  on  differential  pseudolite  navigation. 

1-4  Scope 

The  development  of  an  extended  Kalman  filter  to  produce  the  floating  point  es¬ 
timates  of  carrier-phase  ambiguities  and  the  ambiguity  resolution  techniques  are  the 
focus  of  this  research.  A  simulation  is  used  to  evaluate  the  algorithm’s  integer  ambi¬ 
guity  resolution  performance.  The  scope  of  this  thesis  included  the  simulation  and 
development  of  the  pseudolite  network,  along  with  the  generation  of  error-corrupted 
ranges  between  the  pseudolites  and  the  receivers. 

All  software  development  was  developed  in  the  Matlab  6.5  environment.  An 
evaluation  of  single  versus  widelane  frequency  measurements  is  conducted.  This  the¬ 
sis  includes  the  sensitivity  analysis  of  each  pseudolite  signal  error  source,  including 
the  impact  on  the  ambiguity  resolution  process,  ft  also  explores  new  methodologies 
applied  to  reduce  the  effect  of  errors  and  improve  the  accuracy  of  estimated  ambi¬ 
guities.  Two  of  these  methodologies  were  developed  as  part  of  this  research.  This 
thesis  does  not  address  the  design  and  implementation  of  pseudolite  transmitters 
and  receivers. 
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1.5  Assumptions 

The  following  assumptions  are  made  in  this  thesis: 

a)  Real-time  ambiguity  resolution  is  not  required  because  the  focus  is  the 
augmentation  of  a  post-processed  flight  reference  system. 

b)  The  extended  Kalman  filter  is  not  dependent  upon  a  specific  pseudolite 
system  implementation.  The  pseudolite  can  utilize  transmitters  and  receivers  that 
operate  at  the  GPS  LI  and  L2  frequency  or  in  another  band  such  as  the  S-band.  GPS 
jamming  is  mitigated  for  L-band  transmitters  and  receivers  by  the  use  of  frequency 
translators. 

c)  All  calculations  use  the  Earth-Centered  Earth-Fixed  (ECEF)  frame  and 
World  Geodetic  Systems  1984  (WGS-84)  coordinates. 

d)  No  jamming  analysis  is  required,  because  the  pseudolites  in  a  fielded  system 
would  operate  at  another  frequency. 

e)  The  sources  of  error  present  in  the  code  and  phase  pseudolite  measurements 
are  assumed  to  be  of  similar  characteristics  to  those  available  with  high-quality  GPS 
receivers. 

f)  Carrier-phase  cycle  slip  detection  is  already  accomplished. 

1.6  Thesis  Overview 

Chapter  2  presents  the  background  theory  for  this  research  through  an  in- 
depth  review  of  Kalman  filter  methods,  GPS  fundamentals,  carrier-phase  ambiguity 
resolution,  and  pseudolite  navigation.  The  section  on  Kalman  filtering  includes 
the  derivation  of  an  extended  Kalman  filter,  optimal  smoothing  techniques,  and 
conditional  moment  estimators.  Chapter  3  details  the  error  truth  model,  pseudolite 
filter  models,  and  the  carrier-phase  ambiguity  resolution  techniques  used  in  this 
thesis.  Chapter  4  describes  the  single  run  and  Monte  Carlo  analysis  of  the  effects  of 
each  source  of  measurement  error,  along  with  the  ability  of  filter  enhancements  to 


1-6 


improve  the  accuracy  of  the  position  and  ambiguity  solutions.  Chapter  5  summarizes 
the  results  and  provides  recommendations  for  future  research  on  a  pseudolite-based 
reference  system. 
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II.  Background 


2. 1  Overview 

This  chapter  begins  by  providing  a  basic  overview  of  Kalman  filter  theory, 
including  the  extension  to  extended  Kalman  filter  applications.  The  next  section 
is  on  GPS  operation  and  DGPS  techniques.  From  there  a  section  on  pseudolites 
describes  the  challenges  and  issues  of  pseudolite  navigation.  The  last  section  provides 
the  theory  behind  the  carrier-phase  ambiguity  resolution  techniques  used  in  this 
research. 

2.2  Kalman  Filters 

Deterministic  analysis  has  been  successfully  applied  to  many  systems,  but  is 
not  totally  sufficient  when  applied  to  particular  problems  of  interest.  The  linear 
Kalman  filter  is  an  optimal  recursive  data  processing  algorithm  [19]  that  is  a  common 
tool  technique  that  can  be  applied  when  deterministic  analysis  is  not  sufficient.  The 
optimality  is  based  on  the  assumptions  that  form  the  basis  for  Kalman  filter,  namely, 
an  adequate  model  of  the  real-world  application  in  the  form  of  a  linear  dynamics 
model  driven  by  white  Gaussian  noise  of  known  statistics,  from  which  are  taken 
linear  measurements,  corrupted  by  white  Gaussian  noises  of  known  statistics  [19]. 
The  Kalman  filter  can  produce  very  sub-optimal  results  if  either  the  dynamics  or 
measurement  model  is  an  inadequate  model  of  the  real  world  [34] .  The  Kalman  filter 
is  also  optimal  because  it  incorporates  all  available  measurements,  regardless  of  their 
accuracy,  to  compute  the  estimates  of  the  variables  of  interest  based  on  the  system 
dynamics  and  measurement  models,  the  statistical  description  of  the  system  noises, 
measurement  errors,  and  the  model  uncertainties  [15,  36,  19]. 

When  discrete-time  measurements  are  available,  a  Kalman  filter  includes  both 
a  time  propagation  cycle  and  an  measurement  update  cycle.  The  propagation  cycle 
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computes  an  estimate  of  the  system  state  based  on  its  previous  system  state  and 
its  (imperfect)  dynamics  model.  The  update  cycle  then  uses  the  noise-corrupted 
measurements  to  refine  the  system  state  estimates.  A  complete  derivation  can  be 
found  in  [19]. 

2.2.1  State  and  Measurement  Model  Equations.  The  following  development 
is  similar  to  those  in  references  [15,  19].  The  the  system  dynamics  are  assumed  to 
be  modeled  as  a  linear  system  with  a  state  equation  of  the  form 

x(t)  =  F(t)x(f)  +  B(f)u(f)  +  G(f)w(f)  (2.1) 

where 


x(f)  =  the  n-dimensional  system  state  vector 
F (t)  =  the  n-by-n  state  dynamics  matrix 
B(f)  =  the  n-by-r  control  input  matrix 
u  (f)  =  the  r-dimensional  control  input 
G(f)  =  the  n-by-s  noise  input  matrix 
w(f)  =  the  s-dimensional  dynamics  driving  noise  vector 


Upper  case  bold  letters  indicate  matrices,  lower  case  bold  letters  indicate  vectors,  and 
normal  or  italics  represent  scalar  variables.  Random  vectors  are  denoted  by  boldface 
sans  serif  type.  For  the  purposes  of  this  research  there  are  no  control  inputs,  so  the 
B  and  u  terms  will  be  dropped  from  any  subsequent  equations. 

At  discrete  times  the  solution  to  equation  (2.1)  can  be  written  as: 


x(U+i)  =  $(U+i,U)x(U)  + 


rU+i 

/  $(U+i,r)G(r)d/3(r) 

Ju 


(2.2) 
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where  (3  is  a  vector  valued  Brownian  motion  process  of  diffusion  Q  (f)  [19],  and 
is  the  state  transition  matrix  from  tj  to  t,l+\  and  is  given  by 

=  <I>(A t)  =  eFAt  where  At  =  ti+1  —  U  (2.3) 

which  assumes  a  time-invariant  F  matrix.  The  equivalent  discrete-time  model  is 
expressed  by  the  stochastic  difference  equation  as 


x(fi+1)  =  +  W  d(ti) 


(2.4) 


where 


ti+1 


w  d(U)  =  /  $(ti+i,r)G(r)/3(r) 


(2.5) 


The  discrete-time  white  Gaussian  dynamics  driving  noise  has  the  statistics: 


£{«■((<()}  =  o 


(2,6) 


l'ti+ 1 


E{wd(ti)w^(ti)}  =  Qd(ti)  =  /  +(ti+i,r)G(r)Q(r)Gi  (r)+J  (ti+i,r)dr  (2.7) 


S{wj(ii)wJ(«j)}  -  a.  r.  r  tj 


(2,8) 


Typical  problems  of  interest  are  defined  by  a  continuous-time  dynamics  process  with 
discrete-time  measurements  produced  by  sensors.  Assume  the  measurement  model 
can  be  given  as  a  linear,  discrete-time  system  of  the  form 


z  (U)  =  H  (ii)x(ti)  +  v(ti) 


(2.9) 
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The  statistics  of  the  measurement  corruption  noise  are  described  by 


E{v(U)}  =  0 


R (ti)  for  tj  =  tj 
0  for  ti  7^  tj 


(2.10) 

(2.11) 


The  dynamics  driving  noise  w^(tj)  and  the  measurement  corruption  noise  v(fj)  are 
assumed  to  be  independent,  so 


E{vjd(ti)vT (tj)}  =  0  for  all  ti  and  tj  (2-12) 


2.2.2  Kalman  Filter  Equations.  The  Kalman  filter  propagates  forward 
in  time  from  tf_ 1  to  t~ ,  starting  from  the  last  update  cycle  state  and  covariance 
estimates.  The  superscripts  ”+”  and  ”  denote  the  time  after  a  measurement 
update  and  before  a  measurement  update  respectively.  Propagating  the  filter  from 
ti  to  ti+i  is  equivalent  to  propagating  from  t^i  to  ti.  The  initial  conditions  x(t0) 
and  P(t„)  are  used  in  the  first  propagation  cycle.  The  propagation  cycle  is  given  by 


x(ti)  =  ®(ti,ti-i)x.(tt_1)  (2.13) 

P (ti  )  =  tj_i)P(tJ.i)$T(tj)tj_i)  +  Gd(ti_i)Qd(tj_i)Gj(tj_i)  (2-14) 

When  measurements  are  available,  the  Kalman  filter  is  updated  by: 

A  (U)  =  H(ti)P(tr)HT(ti)  +  R  (U)  (2.15) 

K  (U)  =  P(^r)HT(ti)A(tj)_1  (2.16) 

r (U)  =  z i  -  H (tj)x(tr )  (2.17) 

x(^+)  =  x(t“)  +  K(ti)r(ti)  (2.18) 

P(tf)  =  P (t~)  -  K(A)H(t?:)P (t~)  (2.19) 
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A  properly  designed  filter  has  a  zero-mean  residual  vector  r (ti)  with  the  hlter- 
computed  covariance  A  (ti)  [15].  The  outputs  of  the  Kalman  filter  update  cycle 
are  x(£+)  and  P(tf),  which  are  then  used  in  the  next  propagation  cycle. 

2.3  Extended  Kalman  Filters 

The  linear  Kalman  filter  cannot  be  used  when  either  the  state  dynamics  or 
measurement  model  contains  nonlinearities.  Methods  that  choose  to  ignore  old  data 
due  to  cumulative  errors,  or  that  decrease  the  filter’s  confidence  in  the  adequacy  of 
the  filter  model  have  been  used  to  address  the  problem  of  nonlinearities.  A  better 
way  to  deal  with  nonlinearities  is  through  a  linearization  of  the  measurement  or 
dynamics  model,  thus  enabling  linearized  estimation  techniques. 

A  linearized  Kalman  filter  consists  of  a  first  order  Taylor  series  approximation 
to  the  nonlinear  models,  linearizing  about  a  nominal  trajectory  that  is  normally  pre¬ 
computed.  The  extended  Kalman  filter  (EKF)  differs  from  the  linearized  Kalman 
filter  in  that  it  relinearizes  about  each  state  estimate  as  it  progresses,  enabling  it  to 
handle  larger  degrees  of  nonlinearities  more  adequately.  A  complete  derivation  of 
extended  Kalman  filters  can  be  found  in  reference  [20]. 

2.3.1  State  and  Measurement  Model  Equations.  Following  the  Kalman 
filter  development  in  references  [15,  20],  a  nonlinear  system  dynamics  model  takes 
the  form 

x(f)  =  f[x(f),f]  +  G(t)w(t)  (2.20) 

The  state  dynamics  vector  is  now  defined  to  be  a  possibly  nonlinear  function  of  the 
n-dimensional  state  vector  x(f),  and  of  the  continuous  time,  t,  itself.  The  definitions 
of  the  n-dimensional  state  dynamics  vector  x(i)  and  the  n-by-s  noise  distribution 
matrix  G (t)  are  unchanged  from  those  seen  in  association  with  Equation  2.1.  The 
dynamics  driving  noise  vector  w (t)  is  also  unchanged  and  assumed  to  be  white 
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Gaussian  noise  process  with  mean  and  covariances  defined  by: 


E{w(t)}  =  0 

E{w(t)wJ  (t  +  r)}  =  Q (t)S(r) 


(2.21) 

(2.22) 


where  r  has  units  of  time. 

The  nonlinear  discrete-time  measurement  equation  takes  the  form 


z  (U)  =  h[x(ti),ti]  +  v(ti)  (2.23) 

The  m-dimensional  measurement  vector  z (£j)  is  a  linear  or  nonlinear  function  of  the 
state  vector  and  time  (h[x(fj),  £j]),  corrupted  by  the  linearly  additive  m-dimensional 
discrete-time  noise  input  vector  v(tj).  The  discrete-time  noise  vector  is  unchanged 
from  that  of  the  linear  Kalman  filter. 

2.3.2  State  and  Measurement  Model  Linearization.  If  either  the  system  or 
measurement  model  equations  2.20  and  2.23  are  nonlinear,  they  must  be  linearized 
in  order  to  produce  an  optimal  state  estimate,  to  first  order.  Reference  [20]  uses 
a  perturbation  technique  of  the  state  about  a  nominal  or  reference  trajectory.  The 
dynamics  model  for  this  research  is  linear,  but  the  linearization  of  both  the  dynamics 
model  and  measurement  model  is  presented  for  completeness. 

The  nominal  state  trajectory  can  be  generated  from  the  initial  condition  xn(fo)  = 
x„0  and  the  differential  equation 


xn(t)  =  f  [xn  (f ) ,  f] 


(2.24) 
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which  differs  from  the  nonlinear  state  equation  by  being  deterministic.  The  nominal 
measurements  can  be  defined  in  a  similar  fashion  by 


z  n(ti)  =  h  [xn(t*),ti]  (2.25) 

which  is  also  deterministic.  The  perturbation  state  derivative  Sx(t),  is  formed  by  the 
subtraction  of  the  nominal  trajectory  (2.24)  from  the  system  model  (2.20)  to  give 

6x(t)  =  [x(t)  -  xn(t)}  =  f [x(t),  t]  -  f  [xn(t),  t]  +  G(f)w(f)  (2.26) 

A  Taylor  series  expansion  of  f[x(f),f]  about  xQ(t)  yields 

[x(t)  —  xn(f)]  +  h.o.t.  (2.27) 

x=x„  (t) 

where  ” h.o.t.”  represents  the  higher  order  terms  in  powers  of  [x(t)  —  xn(t)]  greater 
than  one  [20].  When  Equation  2.27  is  substituted  into  Equation  2.26,  the  f[xn(£),i] 
term  is  cancelled  to  produce  the  perturbation  equation.  The  first  order  approxima¬ 
tion  ignores  the  higher  order  terms  which  yields 

Sx(t)  =  F[t;  xn(t)]8x(t)  +  G(t)vj(t)  (2.28) 

This  linearized  dynamics  equation  can  be  implemented  in  a  linearized  Kalman  fil¬ 
ter  with  the  n-by-n  partial  derivative  matrix  F[f;xn(f)]  evaluated  along  a  nominal 
trajectory  and  defined  as 

(2.29) 

X=x„(t) 

This  approximation  is  valid  as  long  as  the  higher  order  terms  of  the  Taylor  series 
expansion  are  negligible. 


F[t;: 


.(*)]= 


dx 


f[x(t),t]  =  f[xn,t]  + 


df{x(t),t] 


dx 
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The  development  of  the  measurement  perturbation  equation  is  formed  in  a 
similar  way.  The  measurement  perturbation  Sz(t)  is  formed  by  the  subtraction  of 
the  nominal  measurement  Equation  2.25  from  the  measurement  model  Equation  2.23 
to  give 

Sz(ti)  =  [z(U)  -  z n(ti)]  =  h[x(ti),ti]  -  h[xn(fj),  ti]  +  v(ti)  (2.30) 

A  Taylor  series  expansion  of  h[x(f;),fj]  about  xn(f)  yields 


h[x(ti),ti 


h[xn(fj),  fj]  + 


dh[x,  U 

dx 


[x(tj)  -  xn(ti)]  +  h.o.t. 

X=Xn(ti) 


(2.31) 


When  Equation  2.31  is  substituted  into  Equation  2.30,  the  h[xn (£*),£*]  term  is  can¬ 
celled  to  produce  the  perturbation  equation.  The  first  order  approximation  ignores 
the  higher  order  terms  which  yields 


5z(ti)  =  H[fj;xn(fj)]5x(^)  +v(^) 


(2.32) 


This  linearized  measurement  equation  can  be  implemented  in  the  linearized  Kalman 
filter  with  the  rn-by-n  partial  derivative  matrix  H[ijjxn(£j)]  evaluated  along  a  nom¬ 
inal  trajectory  and  defined  as: 

(2.33) 

x=xn(ti) 

This  approximation  is  valid  as  long  as  the  higher  order  terms  of  the  Taylor  series 
expansion  in  Equation  2.31  are  negligible.  The  state  and  measurement  perturbation 
equations  are  error  state  representations  which  must  be  added  to  the  nominal  state 
values  to  produce  the  total  state  estimate. 

The  equations  developed  in  this  section  define  the  linearized  Kalman  filter. 
Real-world  measurements  z(U)  are  differenced  with  z n(ti)  computed  via  Equation 
2.25,  and  then  fed  into  a  linear  Kalman  filter  based  on  Equations  2.28  and  2.32,  to 
generate  estimates  of  Sx.(t).  These  can  be  added  to  xn(f),  generated  as  solutions  to 


H[ti;xn(f?:);  = 


<9h[x,  ti 


<9x 
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Equation  2.24,  to  estimate  the  total  states.  It  is  important  to  point  out  that  the 
EKF  relinearizes  the  model  about  the  new  estimate  (x(t+))  and  the  corresponding 
trajectory.  The  relinearization  process  helps  to  validate  the  assumption  that  the 
deviations  from  the  nominal  trajectory  are  sufficiently  small. 

2.3.3  Extended  Kalman  Filter  Equations.  The  extended  Kalman  filter 
propagates  forward  in  time  fj_i  to  f*  by  integrating  from  the  last  update  cycle,  state 
and  covariance  estimates.  The  initial  conditions  x(f 0)  and  P(t0)  are  used  in  the  first 
propagation  cycle.  The  EKF  propagation  equations  are  defined  by: 


x  =  f[x(f|f;_i),f] 


(2.34) 


P(i|*i-r)  =  F[t;x(t|t,_1)]P(t|t,_1)+P(t|ti_1)FT[t;x(t|t,_1)]+G(t)Q(t)GT(t)  (2.35) 

with  t\ti-i  denoting  the  value  of  a  given  variable  at  time  t,  conditioned  on  all  the 
measurements  up  to  and  including  time  fj_i.  The  term  F [f ;  5c(t |ti_1)]  is  the  n-by-n 
partial  derivative  matrix: 


F[t;  x(f|fj_i)' 


<9f[x,u(t),f] 


<9x 


x=x(i|tj_i) 


The  differential  equation  initial  conditions  are  given  by 


(2.36) 


x(ti_ i|ti_ i)  =  x(f+i)  (2.37) 

=P(tti)  (2-38) 


After  integrating  equations  (2.34)  and  (2.35)  to  the  next  sample  time,  the  state  and 
covariance  estimates  are  defined  as: 


X(tj  )  =  X(ti|tj_i) 

P(^i  )  =  P (ti\ti—l) 
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(2.39) 

(2.40) 


The  EKF  incorporates  the  measurements  in  the  following  update  equations: 

K(U)  =  P(t-)HT[ti;  x(fr)]{Hfc;  x(f7)]P(f-)HTk;  x(t“)]  +  R(b)}”1  (2.41) 

x(t+)  =  x(r)  +  K(ti){z i  -  h[x(fr),^]}  (2.42) 

P (t+)  =  P(r)  -  K(ti)H[ti;  x(tf)]P(r)  (2.43) 

2-4  Optimal  Smoothers 

The  traditional  Kalman  filters  runs  forward  in  time — that  is,  it  incorporates 
all  information  available  up  to  and  including  the  current  time.  An  estimator  that 
uses  future  data  to  improve  the  state  and  covariance  estimate  at  the  current  time  is 
called  an  optimal  smoother.  The  three  class  of  smoothers  include  the  fixed- interval, 
fixed-point  and  fixed-lag  smoothers  [20].  The  fixed-interval  smoother  is  the  type  that 
was  used  in  this  research. 

A  fixed-interval  smoother  can  be  conceptualized  by  the  combination  of  two 
filters.  The  first  is  a  traditional  forward-running  Kalman  filter,  and  the  second  is  a 
backward- running  Kalman  filter  that  is  of  inverse  covariance  formulation  [20].  The 
smoothed  estimate  is  formed  by  combining  the  forward  and  backward  state  estimates, 
using  a  weighting  based  on  their  respective  covariance  matrices.  An  equivalent  ap¬ 
proach  to  the  forward-backward  configuration  was  developed  by  Meditch  [20].  The 
filter  requires  that  the  state  and  covariance  matrices  be  stored  both  before  and  after 
measurement  updates.  The  state  transition  matrices  are  also  required.  Once  the 
forward  filter  is  run  through  all  data  until  the  final  time,  the  smoothed  state  esti¬ 
mate  at  the  last  time  increment  is  set  equal  to  the  state  estimate  from  the  forward 
running  filter  after  the  last  measurement  update,  which  is  denoted  by 

*(*/!*/)  =*(*/)  (2-44) 
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Then,  starting  with  the  second-to-last  time  increment,  and  running  backward  in 
time,  the  smoothed  estimate  is  defined  as 

=  x(t+)  +  A(ti)[x(ti+1\tf)  -  x(t"+1)]  (2.45) 

where  the  matrix  A (ti)  is  defined  as  the  smoothing  estimator  gain  matrix  [20],  given 
by 

A(t,)  =  P((+)#T(«i+1,(i)p-1(«-+1)  (2.46) 

In  a  similar  manner,  the  covariance  at  the  final  time  increment  is  defined  as 

P(tf\tf)  =  P(tp  (2.47) 

and  the  covariance  at  every  other  time  increment  is  calculated  by 

P(«il*/)  =  P (V)  +  A((.)[P(«.+1|(/)  -  P(«+)]AT((.)  (2.48) 

Smoothers  outperform  standard  Kalman  filters  particularly  when  the  dynamics  model 
includes  relatively  large  dynamic  driving  noises.  The  more  uncertainty  in  the  model, 
the  greater  the  benefit  from  incorporating  future  information.  At  the  last  time  epoch, 
the  smoothed  estimate  is  equal  to  the  forward  filter  state  estimate.  A  more  rigorous 
derivation  of  smoothers  can  be  found  in  reference  [20]. 

2.5  Global  Positioning  System 

The  Global  Positioning  System  uses  a  constellation  of  medium  earth  orbit 
satellites  to  provide  a  continuous  ranging  source.  The  user  can  calculate  position, 
velocity,  and  time  from  the  received  signal.  Differential  GPS  (DGPS)  is  a  term 
that  includes  many  different  methods  and  techniques  that  result  in  a  greater  accu¬ 
racy  than  stand  alone  GPS.  A  detailed  overview  of  GPS  concepts  can  be  found  in 
references  [21,  26]. 
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2.5.1  GPS  Signal.  The  GPS  signal  contains  both  a  code  and  carrier-phase 
component.  The  GPS  code  that  is  available  to  civilian  users  is  the  Coarse/ Acquisition 
(C/A)  code  while  the  military  also  has  the  precision  (P)  code  (which  is  called  P(Y) 
after  encryption).  The  carrier-phase  frequencies  are  currently  set  at  1575.42  MHz 
and  1227.6  MHz,  which  are  called  the  LI  and  L2  frequencies  respectively  [21].  The 
P(Y)  is  transmitted  on  both  LI  and  L2  while  the  C/A  code  is  only  available  on  the 
LI  frequency.  The  1023  bit  sequence  C/A  code  repeats  every  millisecond  and  the 
P(Y)  repeats  every  7  days  per  satellite.  The  chipping  rates  for  the  C/A  and  P  codes 
are  1.023  MHz  and  10.23  MHz  respectively.  The  code  component  of  the  GPS  signal 
contain  a  pseudorandom  noise  (PRN)  code  that  is  unique  to  each  satellite. 

Civilian  receivers  can  only  track  the  C/A  code  on  the  LI  frequency.  Military 
receivers  that  track  both  the  P(Y)  codes  on  the  LI  and  L2  frequencies  are  called 
dual-frequency  receivers.  Some  civilian  receivers  use  semi-codeless  techniques  that 
can  be  used  to  get  range  information  from  the  P(Y)  code  without  actually  decrypting 
it  [15,  34],  These  high-precision  civilian  receivers  are  used  in  CHAPS  to  get  the  L2 
carrier-phase  information  without  really  tracking  the  P(Y)  code. 

2.5.2  GPS  Measurements.  Typically,  there  are  three  measurements  from 
a  GPS  receiver — code,  doppler,  and  carrier-phase.  The  code  measurement  is  often 
called  a  ”  pseudorange”  because  it  is  the  actual  range  corrupted  by  measurement 
errors  (primarily  the  clock  error).  The  Doppler  measurement  describes  the  frequency 
shift  in  the  signal  due  to  vehicle  (and  clock)  dynamics,  and  the  carrier-phase  can 
be  thought  of  as  an  integrated  Doppler.  The  term  ’’raw”  is  included  to  distinguish 
these  measurements  from  the  navigation  processor  outputs  such  as  position,  velocity, 
and  acceleration.  DGPS  techniques  will  be  distinguished  based  on  the  use  of  code, 
carrier-phase,  or  both. 

2.5.3  Code  Measurements.  The  code  pseudorange  is  true  range  between 
the  satellite  and  user  plus  the  impact  of  a  number  of  error  sources.  The  pseudo- 
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range  is  calculated  as  the  time  difference  between  the  transmission  and  reception 
time  multiplied  by  speed  of  light  (to  give  the  range  in  meters).  The  pseudorange 
measurement  can  be  expressed  as: 

p  —  r  +  c(Stu  -  Stsv)  +T  +  I  +  mp  +  vp  (2.49) 


where 

p  =  GPS  pseudorange  measurement  (meters) 
r  =  true  range  from  the  user  to  satellite  (meters) 
c  =  speed  of  light  (meters  /  second) 

Stu  =  receiver  (user)  clock  error  (seconds) 

8tsv  =  transmitter  (satellite  vehicle)  clock  error  (seconds) 

T  =  errors  due  to  tropospheric  delay  (meters) 

/  =  errors  due  to  ionospheric  delay  (meters) 
mp  =  errors  due  to  pseudorange  multipath  (meters) 
vp  =  errors  in  pseudorange  clue  to  receiver  noise  (meters) 

2.5.4  Carrier- Phase  Measurements.  The  carrier-phase  of  the  received  sig¬ 
nal  can  also  be  used  for  positioning,  especially  when  high  precision  is  required.  The 
carrier-phase  measurement  can  be  expressed  in  cycles  as: 

0  =  A-1(r  +  c(Stu  -  Stsv )  +  T  -  I  +  rricj,  +  v $)  +  N  (2.50) 
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where 


(/)=  carrier-phase  measurement  (cycles) 

A=  carrier-phase  wavelength  (meters  /  cycle) 

N—  carrier-phase  integer  ambiguity  (cycles) 

T—  errors  clue  to  tropospheric  delay  (meters) 

1=  errors  due  to  ionospheric  delay  (meters) 

77^=  errors  due  to  carrier-phase  measurement  multipath  (meters) 

errors  in  carrier-phase  measurement  due  to  receiver  noise  (meters) 

The  rest  of  the  terms  were  previously  defined  in  Equation  (2.49)  except  that  the 
measurement  noise  and  multipath  are  different  and  significantly  less  for  the  carrier- 
phase  than  for  the  code.  Some  sources  of  error  do  not  affect  the  code  and  carrier- 
phase  measurement  in  the  same  manner.  The  sign  on  the  ionospheric  delay  term 
is  opposite  from  the  code  measurement  equation,  because  the  ionosphere  advances 
a  carrier-phase  measurement,  but  delays  a  code  measurement.  This  phenomenon  is 
called  code-carrier  divergence.  Tropospheric  delay  affects  both  the  code  and  carrier- 
phase  by  the  same  magnitude. 

The  carrier-phase  integer  ambiguity  term  is  an  error  source  that  is  present  in 
carrier-phase  measurements,  but  not  in  code  measurements.  The  ambiguity  term 
represents  the  unknown  number  of  carrier-cycles  present  at  the  start  of  the  signal 
integration  [28].  The  carrier  signal  can  be  broken  up  into  two  segments  that  are 
separated  by  the  point  in  time  that  signal  integration  starts.  The  first  segment 
consists  of  the  unknown  integer  number  of  cycles  up  to  the  point  of  signal  integration. 
The  second  segment  consists  of  the  measured  carrier-phase  which  is  not  constrained 
to  be  an  integer.  A  high  level  of  accuracy  is  attained  when  the  unknown  number  of 
cycles  before  signal  integration  is  determined. 
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2.5.5  Single  Differencing.  Differential  GPS  uses  linear  combinations  of 
observations  (code  or  carrier  measurements)  between  receivers,  satellites,  or  times 
to  reduce  the  effect  of  some  errors  [22],  A  single  difference  can  be  between  two 
satellites  (V)  or  between  two  receivers  (A).  Figure  2.1  depicts  the  concept  of  a 
single  difference  between  two  receivers. 


/  \ 
/  \ 


A 


B 


Figure  2.1  Single  Difference  GPS  Between  Receivers  A  and  B  and  Satellite  k 


The  single-differenced  carrier-phase  measurement  between  two  receivers  corre¬ 
sponding  to  the  above  figure  is  defined  as 

A  4ab  =  <Pb  (2.51) 

where  <pkA  is  the  phase  measurement  between  receiver  A  and  satellite  /c,  and  (j)kB  is 
the  phase  measurement  between  receiver  B  and  satellite  k. 

This  type  of  difference  eliminates  the  satellite  clock  error  and  reduces  the 
atmospheric  errors.  Combining  the  carrier-phase  measurement  Equation  2.50  with 
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the  single  difference  Equation  2.51  yields 


A  <&B  =  A-‘[  r\  +  c  «,  -  UtA)  +  TLA-l\  +  m\A  +  v*A]  +  N\ 

-  A-‘[r|  +  c«„  -  6&J  +T%- 4  +  m\B  +  «JB]  +  JV*  (2.52) 

Combining  like  terms  yields 

a  <&b  =  A-'[(q  -  4)  +  ■=«,  -  <,)  -  +  (d  -  d) 

-  (d  -  4)  +  «4  -  ™*b)  +  «.  -  »*b)1  +  (d  -  JVs)  (2-53) 

After  the  satellite  clock  term  is  eliminated  (because  the  measurements  are  syn¬ 
chronous  and  the  satellite  clock  error  is  the  same  for  both),  differences  are  represented 
as  (A),  and  the  above  equation  can  be  written  as 

A 4>kAB  =  A-X(A rkAB  +  cA5tkUAB  +  A T\b  -  l\B  +  A m\AB  +  vkAB)  +  A NkAB  (2.54) 

The  integer  value  ANab  is  the  difference  in  the  carrier-phase  ambiguity  between  the 
two  receivers’  measurements. 

2.5.6  Double  Differencing.  A  double  difference  is  the  combination  of  a 
single  difference  between  satellites  (transmitters)  and  a  single  difference  between 
receivers.  Because  a  single  difference  between  receivers  cancels  the  satellite  clock 
error  and  a  single  difference  between  satellites  cancels  the  receiver  clock  error,  the 
double  difference  cancels  both  clock  error  terms.  Figure  2.2  displays  the  concept  of 
a  double  difference. 
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Figure  2.2  Double  Difference  Between  Satellites  j  and  k  with  Receivers  A  and  B 

The  phase  measurement  will  be  used  in  the  following  example.  The  double 
differenced  carrier-phase  measurement  is  defined  as 

=  A  (f)kAB  —  A  ^AB  (2.55) 

After  the  single  differenced  phase  Equation  2.54  is  substituted  in  the  above  equation, 
it  yields 

A V<t>kABj  =  A-'(A rkAB  +  cA5tkUAB  +  A TkAB  -  A l\B  +  A m\B  +  AvkAB  +  A NkB 

~  [A-1  (A rJAB  +  cA5tJUAB  +  A T\b  -  A I\B  +  A m\AB  +  Av\ab  +  A NJAB\ 

(2.56) 

When  the  user  clock  error  term  is  cancelled  and  the  double  difference  operator  (AV) 
is  used  to  express  the  double  difference  error  terms,  the  equation  can  be  written  as 

A^Jb  =  X-\AVr%+AVT%-AVI%+AVm%B+AVv%B)+AVN%  (2.57) 
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Differencing  reduces  the  effects  of  frequency  correlated  errors  (such  as  the  atmo¬ 
spheric  errors)  at  the  expense  of  increasing  the  effects  of  uncorrelated  frequency 
errors  (such  as  measurement  noise  and  multipath).  The  single  difference  increases 
the  magnitude  of  the  noise  and  multipath  by  a  factor  of  (\/2)  and  the  double  differ¬ 
ence  increases  the  magnitude  by  a  factor  of  2.  Although  the  integer  ambiguity  term 
(AViV$i)  is  different  from  the  ambiguity  term  from  the  observation  equation,  it  has 
maintained  its  integer  nature. 

The  double  differenced  code  measurement  can  be  adapted  from  Equation  2.57 
by  dropping  the  ambiguity  terms  and  expressing  the  range  in  terms  of  meters,  yield¬ 
ing 


AV&  =  AVr*,  +  AVrg,  -  A VI%  +  AVm«  B  +  AV»«  B  (2.58) 

It  is  important  to  note  that  double-difference  code  measurements  are  not  typ¬ 
ically  used.  Rather,  single  difference  measurements  are  used  and  the  receiver  clock 
error  is  estimated  directly. 

2.5.1  Widelane  Measurements.  When  two  GPS  frequencies  are  available, 
linear  combinations  can  be  formed  to  create  new  virtual  frequencies.  The  widelane 
measurement  is  commonly  used  in  DGPS  applications  and  can  be  written  as  [28] 

4>wl  =  4>li  —  4>L2  (2.59) 

The  widelane  measurement  has  a  wavelength  of  approximately  86.19  cm,  while  the 
wavelengths  of  signals  at  the  LI  and  L2  frequency  are  19.03  cm  and  24.42  cm, 
respectively.  Although  the  integer  ambiguity  value  for  a  widelane  measurement  is 
different  from  either  of  the  values  for  the  LI  or  L2  frequency,  it  is  still  an  integer. 

Table  2.1  shows  a  comparison  between  error  sources  for  a  widelane  versus  single 
frequency  in  terms  of  cycles,  where  c  represents  the  speed  of  light  and  fi  and  /2  are 
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Table  2.1  Comparison  between  LI  and  wi delane  (WL)  phase  errors  (cycles) 


Error  LI  error  (cycles)  WL  error  (cycles)  WL/L1  ratio 


SV  Position 

^-VA  5psv 

a^VA<W 

Ai  ^ 
^ WL 

Ai  ^ 

V l  „ . 

0.221 

Troposphere 

fVA  T 

Ai 

.  1  VAT 

XWL 

0.221 

Ionosphere 

(/l'/2)VA  / 

c/1/2 

(/1-/2) 

h 

«  0.283 

LI  Multipath 

XfVAm, 

vVA^ 

1 

LI  Noise 

VVAW 

VVA^ 

1 

L2  Multipath 

AVAm0 

VVA^ 

1 

L2  Noise 

VVAW 

vVAW 

1 

the  frequencies  of  LI  and  L2  respectively  [28].  The  widelane  significantly  reduces 
the  impact  of  the  correlated  error  sources  (i.e. ,  the  satellite  and  receiver  position 
errors  and  the  atmospheric  errors).  Measurement  noise  and  multipath  are  frequency 
uncorrelated  and  thus  not  reduced  by  a  widelane  implementation. 

Although  the  widelane  measurement  reduces  some  of  the  error  sources  when 
expressed  in  cycles,  it  actually  increases  some  of  the  error  sources  when  expressed 
in  meters.  To  convert  the  widelane  to  LI  ratio,  a  conversion  of  «  4.529  is 
multiplied  by  the  ratio  give  in  the  table.  This  means  that,  when  the  multipath 
and  measurement  noise  are  not  affected  in  terms  of  cycles,  the  effect  is  amplified  by 
a  factor  of  approximately  4.529  when  expressed  in  meters.  The  satellite  position, 
receiver  position,  and  tropospheric  errors  are  reduced  when  expressed  in  cycles,  but 
are  unaffected  when  expressed  in  terms  of  meters  of  error.  The  ionospheric  error 
is  slightly  increased  in  terms  of  meters  for  a  widelane  measurement,  but  this  error 
source  will  be  ignored  for  pseudolite  applications. 

The  longer  wavelength  of  a  widelane  measurement  reduces  the  number  of  can¬ 
didate  ambiguity  sets  that  are  generated  within  the  ambiguity  resolution  search  space 
(see  Section  2.5.8).  This  makes  ambiguity  resolution  easier  and  more  reliable,  but 
with  a  decreased  accuracy  when  compared  to  a  single  frequency,  because  the  errors 
in  the  widelane  phase  measurements  are  actually  larger  than  the  single  frequency 
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case.  A  more  detailed  discussion  on  linear  combinations  of  measurements  can  be 
found  in  references  [22,  28]. 

2.5.8  Carrier-Phase  Ambiguity  Resolution.  Carrier- phase  ambiguity  reso¬ 
lution  is  the  process  of  selecting  the  correct  integer  value  for  the  phase  ambiguity.  It 
is  not  always  possible  to  perform  ambiguity  resolution,  and  sometimes  the  wrong  in¬ 
teger  is  chosen,  producing  erroneous  results.  Ambiguity  resolution  generally  consists 
of  two  primary  operations  [27].  The  first  is  to  create  the  ambiguity  search  space  by 
the  generation  of  candidate  ambiguity  sets.  The  second  operation  is  the  selection  of 
the  correct  ambiguity  set.  The  next  two  sections  cover  the  algorithms  used  in  this 
research. 


2.5.8. 1  Ambiguity  Set  Generation.  Ambiguity  set  generation  can  be 
characterized  by  being  either  a  position-domain,  measurement-domain,  or  ambiguity- 
centered  approach.  This  section  provides  an  overview  of  the  Least  squares  AMbiguity 
Decorrelation  Adjustment  (LAMBDA)  developed  by  Teunissen  [31]  and  the  Fast 
Ambiguity  Search  Filter  developed  by  Chen  and  Lapachelle  [4] 

LAMBDA  is  not  a  set  generation  technique,  but  rather  a  search  space  trans¬ 
formation  technique.  The  ambiguity  estimates  of  the  floating-point  solution  contain 
a  high  degree  of  correlation,  which  makes  ambiguity  resolution  difficult.  LAMBDA 
reduces  the  correlation  of  the  ambiguity  estimates  to  enable  a  fast  and  efficient 
search  [31].  Teunissen  referred  to  the  ambiguity  transformation  process  as  a  ”Z- 
transformation” ,  not  to  be  confused  with  the  Z  transformation  of  discrete  time  signal 
processing.  The  LAMBDA  method  preserves  the  volume  of  the  search  space  while 
also  maintaining  the  integer  nature  of  the  ambiguities.  The  Z-Transform  is  defined 
as: 

z  =  Z2x  z  =  Z7x  Qz  =  ZtQxZ  (2.60) 
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where 


x,  z  =  untransformed  and  transformed  ambiguities,  respectively 
x,  z  =  untransformed  and  transformed  ambiguity  estimates,  respectively 
Z  =  Z-Transformation  matrix 

Qx,  Qz  =  untransformed  and  transformed  covariance  matrix,  respectively 

The  transformation  technique  can  be  constructed  for  the  simple  two-dimensional 
case,  starting  with  the  untransformed  covariance  ambiguity  covariance  matrix  given 
by 

ai  a  12 

Qx  =  ;  (2.6i) 

°21  a2 

In  this  case,  the  Z-Transformation  matrices  can  be  defined  as 

_  1  int(-ai2a2“2) 

0  1 

or  alternatively 

r 1  o  / 

z  l  =  (2.63) 

int(-cr2icr1  )  1 

where  either  the  upper  (Z\ )  or  lower  (Zj)  diagonal  form  can  be  used.  The  rounding 
of  the  off-diagonal  terms  (denoted  by  ”  int”  )  to  integers  is  necessary  to  maintain  the 
integer  nature  of  the  ambiguities. 

Rizos  and  Han  [29]  developed  an  efficient  method  for  high-order  Z-transformations. 
The  first  step  is  to  perform  an  upper  triangular  factorization  of  the  ambiguity  co- 


(2.62) 
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variance  matrix  by: 


Qx  =  UrD^Uf 


(2.64) 


where  Ui  is  an  upper  triangular  matrix  and  D[r,  is  a  diagonal  matrix. 

In  the  next  step  an  intermediate  transformation  matrix  is  computed  by  first 
rounding  the  elements  of  Ui  to  integer  values  and  then  taking  the  inverse. 

Z  Vl  =  QUiU)-1  (2.65) 

An  intermediate  covariance  matrix  (Q^  )  uses  the  intermediate  transformation  ma¬ 
trix  (ZrJl)  by 

Q  zVl  =  ZuMtZu,  (2-66) 

This  process  is  repeated  again,  except  with  a  lower  triangular  factorization  given  by 

QzLl  =LiDLlLf  (2.67) 

where  Lx  is  a  lower  triangular  matrix  and  D^1  is  a  diagonal  matrix. 

Again  the  matrix  is  first  rounded  and  then  the  inverse  is  taken. 

ZLl  =  ([LiU)"1  (2.68) 

The  intermediate  covariance  is  calculated  by: 

QzLl  =  Zl1Qzc/iZLi  (2.69) 

The  process  is  repeated  until  Equations  2.65  and  2.68  result  in  identity  matrices. 

The  Z-transformation  is  then  given  by 

Z  =  ZLk_1Z[/k_1...ZL2Z[/2ZLlZ;y1  (2.70) 
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Because  the  LAMBDA  method  does  not  actually  generate  an  ambiguity  search  space, 
other  ambiguity  resolution  techniques  are  used  on  the  transformed  estimates  and 
resulting  covariances.  The  Fast  Ambiguity  Search  Filter  (FASF)  is  the  method  used 
in  this  thesis. 

The  FASF  method  was  developed  by  Chen  and  Lachapellc  [4]  as  an  efficient 
approach  for  ambiguity  resolution.  FASF  operates  recursively  and  takes  advantage 
of  the  related  nature  of  ambiguity  ranges.  Along  with  the  ambiguity  state  and 
covariance,  FASF  requires  the  constant  k  to  be  given  to  define  the  search  space  for 
the  ambiguity  VAiYint  by 


x„  -  kan  <  VATVint  <  xn  +  kan  (2.71) 

with  xn  and  crn  representing  the  floating-point  estimate  and  standard  deviation  of 
the  nth  ambiguity.  Conditional  ambiguity  estimates  and  their  associated  covariances 
are  determined  with  the  condition  that  the  first  ambiguity  is  correct.  The  process 
is  recursive  with  each  new  ambiguity  calculated  in  the  same  manner  [4],  The  new 
conditional  state  estimate  and  covariance  are  defined  by: 

x  =  x  -  p n(xn  -  VA Nint)/al  (2.72) 

Pi  =  Px-  (PnPD/^n  (2-73) 


where 

x  =  unconditioned  estimated  parameter  vector 
Pi  =  unconditioned  estimated  parameter  covariance  matrix 
x  =  estimated  parameter  vector  conditioned  on  xn  =  VA N-mt 
P;T  =  covariance  matrix  conditioned  on  xn  =  VA N[nt 
pn  =  nth  column  of  the  covariance  matrix  P £ 

a2  =  scalar  variance  of  the  nth  parameter  (taken  from  the  diagonal  of  Prc) 
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FASF  can  be  performed  in  either  the  normal  or  Z-transformed  domain.  If 
FASF  is  performed  in  the  Z-transform  domain  (as  done  in  this  thesis  research),  it 
must  be  transformed  back  into  the  original  ambiguity  domain. 

2. 5. 8. 2  Ambiguity  Set  Determination.  Ambiguity  set  determination 
is  the  second  operation  in  the  ambiguity  resolution  process.  A  common  technique 
used  in  set  discrimination  is  the  use  of  a  ratio  test  using  a  sum  of  squared  residuals. 
The  correct  set  of  ambiguities  typically  has  the  smallest  residuals.  By  comparing  the 
sum  of  squared  residuals  the  correct  set  generally  stands  out.  This  techniques  can 
be  further  broken  down  into  two  methods  based  on  how  the  residuals  are  formed. 

The  first  method  is  to  define  a  residual  as  the  difference  between  the  floating¬ 
point  ambiguity  state  estimate  and  each  candidate  set  [24],  When  this  approach  is 
used,  the  weighted  sum  of  squared  residuals  is  expressed  by: 

=  (Afloat  ~  (Afloat  ~  Xi)  (2.74) 

where  5q  is  the  ith  candidate  integer  ambiguity  vector  and  0*  is  the  sum  of  squared 
residuals  for  the  ith  candidate  ambiguity  vector.  This  is  the  method  used  in  this 
research. 

The  second  method  involves  defining  residuals  based  on  the  difference  between 
the  measurement  and  measurement  prediction  such  as 

r*  =  Zj  —  h[x(fj);  ti]  (2.75) 

This  requires  that  a  filter  operates  on  each  candidate  ambiguity  set  to  determine  the 
best  fit  for  discrimination.  The  sum  of  squared  residuals  for  this  convention  is  given 
as 

=  rTA_1r  (2.76) 
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where  r  is  a  residual  based  on  equation  (2.75),  and  A  is  defined  to  be  the  measure¬ 
ment  covariance  [24], 


The  ratio  test  can  be  applied  to  either  residual  convention  by  comparing  the 
magnitude  of  the  best  (smallest)  and  second  best  (second  smallest)  sum  of  squared 
residual  terms  by 


ratio  = 


Di(2ndbest) 


(2.77) 


Vti(best ) 

The  ambiguity  set  corresponding  to  the  one  selected  as  best  is  determined  to  be  the 
correct  one  if  the  ratio  is  consistently  large  (typically,  greater  than  2). 


2.6  Pseudolites 

The  term  pseudolite  is  short  for  ’’pseudo-satellite”,  refering  to  ground-based 
GPS-like  transmitters.  Pseudolites  have  the  flexibility  to  vary  the  location,  power, 
and  frequency  of  the  transmitter.  Pseudolites  also  are  able  to  provide  signals  for 
navigation  purposes  in  adverse  environments  such  as  open-pit  mining,  where  GPS 
signals  are  often  block  by  the  steep  sides  of  the  pit.  Many  of  the  assumptions  that 
are  made  with  GPS  navigation  cannot  be  applied  to  pseudolites,  as  will  be  seen. 
This  section  begins  with  a  discussion  of  differences  between  GPS  and  pseudolite 
navigation,  then  presents  typical  pseudolite  applications,  and  ends  with  descriptions 
of  the  problems  and  sources  of  error  in  pseudolite  measurements. 

2. 6. 1  GPS -Pseudolite  Differences.  Many  of  the  assumptions  that  are  used 
in  GPS  processing  cannot  be  made  for  pseudolite  navigation.  These  include: 

•  The  expected  ranges  for  pseudolites  are  much  more  dynamic  than  for  GPS 
operation  and  will  affect  receiver  power  levels 

•  When  a  static  reference  receiver  is  used,  there  is  not  any  relative  motion  be¬ 
tween  the  reference  receiver  and  each  pseudolite  like  there  is  between  a  ref¬ 
erence  receiver  and  the  orbiting  GPS  satellites.  This  results  in  measurement 
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biases  due  to  pseudolite  location  errors  that  do  not  average  out  over  time. 
Also  the  multipath  error  between  pseudolites  and  the  reference  receiver  will 
have  stronger  time  correlations  than  the  multipath  experienced  at  the  mobile 
receiver,  when  in  motion. 

•  Pseudolites  do  not  have  to  operate  at  the  GPS  LI  and  L2  frequencies.  They 
also  do  not  have  to  use  the  same  code  sequence  or  chipping  rates  that  GPS 
satellites  use. 

•  Due  to  the  short  ranges  between  pseudolite  transmitters  and  receivers  the 
measurement  model  is  more  nonlinear  than  for  GPS  operation. 

•  Pseudolites  do  not  have  an  orbital  or  ephemeris  error,  but  rather  a  position 
error  that  is  dependent  on  the  type  of  surveying  accuracy  used  to  estimate  the 
phase  center  of  the  pseudolite  antenna. 

•  A  pseudolite  signal  will  not  travel  through  the  ionosphere,  so  any  error  terms 
associated  with  ionospheric  delay  can  be  ignored. 

•  Depending  on  the  relation  of  the  mobile  receiver  to  the  reference  receiver, 
single  and  double  differencing  will  not  reduce  tropospheric  error  as  much  as 
with  GPS.  This  is  similar  to  extremely  long  baselines  in  GPS  processing. 

2.6.2  Pseudolite  Applications.  The  four  categories  of  pseudolite  applica¬ 
tions  include  direct  positioning,  digital  data  transmission,  carrier-phase  ambiguity 
resolution,  and  as  a  differential  reference  station  [7].  Direct  positioning  using  a 
network  of  pseudolites  is  the  application  addressed  in  this  research. 

Pseudolite  direct  positioning  can  be  accomplished  with  both  the  code  and 
carrier-phase  measurements  in  a  similar  manner  as  for  conventional  GPS  positioning. 
The  majority  of  work  with  pseudolites  has  been  concerned  with  the  augmentation  of 
GPS  or  GPS/INS.  Pseudolites  can  improve  the  overall  geometry  of  the  augmented 
system,  providing  greater  positioning  accuracy,  reliability,  availability,  continuity, 
and  integrity  monitoring  [23].  They  can  also  be  used  as  a  sole  source  signal  of 
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navigation.  Additionally,  GPS  signals  are  typically  weak  or  not  present  indoors,  and 
pseudolites  can  be  used  to  provide  an  indoor  navigation  source. 

Digital  data  transmission  is  also  possible  for  pseudolite  transmitters.  One  of 
the  advantages  of  using  a  pseudolite  for  data  transmission  to  a  GPS  receiver  is  that 
only  minor  changes  are  needed  in  the  receiver  [7].  GPS  reference  data  has  frequently 
been  proposed  for  transmission  via  pseudolites  [10,  30]. 

Pseudolites  can  assist  and  speed  up  the  carrier-phase  ambiguity  resolution  in 
a  GPS  system  augmented  by  pseudolites.  This  is  accomplished  through  the  large 
changes  in  geometries  of  the  pseudolite  signal  [23].  An  example  of  this  is  the  Kine¬ 
matic  GPS  Landing  Systm  (KGLS)  at  Stanford  [7]. 

When  a  pseudolite  rebroadcasts  a  coherent  replica  of  received  GPS  signals,  it 
is  known  as  a  differential  reference  station  [7] .  The  difference  between  the  direct  and 
reflected  signal  can  be  used  for  navigation  purposes. 

2.6.3  Signal  Interference  and  Near-Far  Problem.  One  of  the  largest  issues 
facing  practical  pseudolite  applications  is  the  signal  interference  and  the  associated 
near-far  problem.  While  the  distance  from  a  given  receiver  to  any  GPS  satellite  is 
relatively  constant,  the  ranges  between  a  pseudolite  and  its  receiver  vary  greatly. 
The  large  dynamic  difference  in  ranges  result  in  large  differences  in  received  power 
levels.  This  can  cause  the  automatic  gain  control  in  a  receiver  to  adjust  to  the  highest 
powered  signal,  which  effectively  jams  all  other  pseudolites. 

Pseudolites  have  both  a  ’’near”  and  a  ’’far”  radius  that  define  its  usable  area. 
A  pseudolite  will  jam  all  other  pseudolites  within  its  near  radius.  The  far  radius  is 
the  distance  within  which  a  receiver  must  stay  to  maintain  lock  on  that  pseudolite. 
The  near  and  far  radii  are  a  function  of  transmission  power,  so  increasing  or  decreas¬ 
ing  power  will  increase  or  decrease  the  near  and  far  radii  by  the  same  ratio.  The 
relationship  between  the  near  and  far  radius  is  given  as  a  ratio  that  is  generally  ac¬ 
cepted  to  be  1/10  for  practical  systems  [7],  although  this  can  vary  depending  on  the 
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cross-correlation  of  the  codes.  An  example  of  the  near-far  radii  is  shown  in  Figure 
2.3. 


Figure  2.3  Near-Far  Problem 

The  various  techniques  that  have  been  proposed  to  reduce  the  near-far  prob¬ 
lem  can  be  grouped  into  three  categories — Time  Division  Multiple  Access  (TDMA), 
Frequency  Division  Multiple  Access  (FDMA),  and  Code  Division  Multiple  Access 
(CDMA). 

TDMA  can  be  accomplished  through  the  pulsing  of  the  pseudolites.  If  pulsed 
pseudolites  are  operated  at  greater  than  20-25  percent  of  a  duty  cycle,  then  the  GPS 
signal  will  be  jammed  [5].  It  has  been  proposed  to  operate  two  pseudolites  each  puls¬ 
ing  at  10-12.5  percent  of  the  duty  cycle  to  facilitate  an  integrated  GPS/Pseudolite 
navigation  system  [5].  This  arrangement  still  only  allows  the  use  of  two  pseudolites 
if  the  GPS  signal  is  also  desired.  If  the  GPS  signal  is  not  of  interest,  10  pseudolites 
could  be  used  (given  a  10  percent  duty  cycle  each). 

The  second  technique  for  interference  reduction,  FDMA,  can  be  implemented 
by  modifying  GPS  signals  with  small  frequency  offsets.  Elrod  et.  al  [11]  suggested 
offsetting  the  frequency  to  the  first  null  of  the  GPS  satellite  signal  in  order  to  reduce 
cross-correlation  with  the  GPS  signal.  It  resembles  a  large  doppler  offset  that  most 
receivers  can  handle. 
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CDMA  has  been  demonstrated  through  concatenations  of  C/A  codes.  Ndili  [23] 
showed  through  simulation  that  a  code  length  of  4092  (4  times  that  of  C/A  code) 
would  provide  a  6  dB  enhancement,  and  thus  double  the  far  radius  while  maintaining 
the  same  near  radius.  By  combining  20  C/A  codes  for  a  length  of  20460  in  addition 
to  operating  at  a  P-code  chipping  rate  would  add  a  23  dB  enhancement.  The  longer 
the  code  length  and  higher  the  chipping  rate,  the  larger  the  near-far  ratio. 

2.6.4  Sources  of  Error.  Pseudolite  error  can  be  broken  up  into  measure¬ 
ment  and  measurement  model  errors.  Both  of  these  will  affect  the  accuracy  of  a 
pseudolite  system  because  the  residual  term  is  formed  by  subtracting  the  measure¬ 
ment  prediction  from  the  measurement  as  shown  by  the  following  equation: 

r*  =  Zj  —  h[x(fj);  tf\  (2.78) 

The  next  two  subsections  describe  the  errors  present  in  pseudolite  measurements  and 
measurement  models. 

2.6.4- 1  Pseudolite  Measurement  Errors.  Pseudolite  signal  errors  are 
reduced  less  by  double  differencing  than  for  the  analogous  GPS  equations  due  to  a 
different  geometric  configuration  [6].  The  errors  that  remain  in  a  pseudolite  carrier- 
phase  measurement  after  a  double  difference  operation  are  the  measurement  noise, 
multipath,  and  residual  tropospheric  error  (i.e.,  the  error  after  a  tropospheric  model 
has  been  applied). 

The  measurement  noise  associated  with  pseudolites  is  determined  by  the  qual¬ 
ity  of  the  receiver  (just  like  for  a  GPS  measurement).  Along  with  proper  modelling 
in  the  navigation  filter,  improving  the  receiver  design  is  one  of  the  few  ways  to  reduce 
the  effect  of  measurement  noise. 

Multipath  can  be  considered  a  dominant  error  source  in  pseudolite  applica¬ 
tions  [35].  Even  after  multipath  mitigating  techniques  (such  as  antenna  placement 
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and  chokering  antennas)  are  used  in  a  pseudolite  system,  the  multipath  error  is  gen¬ 
erally  larger  for  a  pseudolite  signal  than  for  a  satellite  signal  [8].  This  is  due  to  the 
relative  geometries  in  the  transmitter-receiver  setup  in  a  pseudolite  network. 

Because,  unlike  satellites,  the  transmitters  do  not  move  relative  to  the  refer¬ 
ence  receiver,  time- invariant  or  standing  multipath  is  a  concern  for  pseudolite  appli¬ 
cations.  This  contributes  to  a  multipath  error  that  is  more  difficult  to  handle  than 
the  multipath  error  associated  with  satellite  signals.  The  navigation  filter  assumes 
that  all  errors  sources  are  white  (uncorrelated  in  time)  and  the  more  time-invariant 
the  multipath  becomes,  the  further  this  assumption  is  broken.  Careful  calibration 
is  required  to  estimate  and  remove  this  constant  error  from  the  corresponding  mea¬ 
surements.  The  use  of  carrier-phase  measurements  and  antenna  design  have  been 
shown  to  help  reduce  multipath  [35].  Multipath  affects  code  measurements  to  a 
higher  degree  than  carrier-phase  measurements  for  both  a  pseudolite  or  satellite 
source.  Antenna  gain  shaping  helps  to  reduce  multipath  by  adjusting  the  gain  in  the 
direction  of  large  reflectors  [18,  35]. 

The  residual  tropospheric  error  that  exists  after  a  tropospheric  model  is  applied 
cannot  be  ignored  in  precise  pseudolite  applications.  The  amount  that  single  and 
double  differencing  reduce  the  effect  of  tropospheric  delay  in  GPS  operation  is  a 
function  of  the  baseline  difference  in  the  mobile  and  reference  receiver  positions. 
Pseudolite  applications  are  equivalent  to  a  very  large  baseline  for  which  differencing 
may  reduce,  but  not  significantly  remove,  tropospheric  delay. 

2. 6. 4-2  Pseudolite  Measurement  Model  Errors.  Pseudolite  measure¬ 
ment  model  errors  include  the  effect  of  position  errors  in  the  location  of  the  pseudo- 
lites  and  reference  receiver  in  addition  to  the  error  clue  to  linear  approximations  in 
the  measurement  model. 

The  source  of  errors  from  the  imprecise  locations  of  the  pseudolite  transmitters 
and  reference  receivers  are  analogous  to  the  ephemeris  or  orbital  errors  in  GPS 
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satellite  locations.  Like  tropospheric  errors,  these  position  errors  are  not  reduced  for 
pseudolites  as  much  as  GPS  for  single  and  double  differencing. 

For  outdoor  pseudolite  applications,  static  surveying  techniques  that  use  carrier- 
phase  DGPS  can  be  used  to  solve  for  the  positions  of  the  pseudolites  to  within 
centimeters.  Because  GPS  is  generally  not  available  indoors,  this  cannot  be  used 
to  estimate  the  locations  of  pseudolite  transmitters  indoors  precisely.  Changdon 
et.al.  [2]  presented  a  method  to  calculate  the  pseudolite  positions  using  only  the 
user’s  position  information  and  the  pseudolite  signals.  This  is  advantageous  because 
the  location  found  was  the  phase  center  of  the  antenna  and  not  just  the  physical 
center. 

One  of  the  biggest  differences  between  GPS  and  pseudolite  navigation  is  the 
typical  ranges  between  transmitters  and  receivers.  GPS  signals  travel  on  the  order 
of  20,000  kilometers  or  more,  while  pseudolite  signal  ranges  could  be  measured  in 
meters  (depending  on  signal  power).  As  the  ranges  in  pseudolite  navigation  become 
shorter,  the  signal  waveform  becomes  more  spherical  than  planar.  Figure  2.4  depicts 
this  relationship  with  a  planar  signal  from  a  GPS  satellite  and  a  spherical  signal 
from  a  pseudolite.  In  reality  a  GPS  signal  is  spherical,  but  the  radius  is  so  large  that 
it  is  essentially  planar  for  a  GPS  user. 


Figure  2.4  Spherical  and  Planar  Wavefronts 
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The  measurement  model  equation  for  the  extended  Kalman  filter  (EKF)  is 
nonlinear  for  both  GPS  and  pseudolite  navigation.  An  EKF  does  a  linearization 
of  the  nonlinear  measurement  equation  by  using  a  first  order  Taylor  series  approx¬ 
imation.  As  the  waveforms  become  more  spherical,  the  measurement  nonlinearity 
becomes  more  severe  and  the  first  order  approximation  becomes  more  inadequate. 
For  GPS  signals,  the  approximation  error  is  small  enough  to  be  ignored,  but  pseu¬ 
dolite  navigation  is  different  and  requires  care  in  handling  the  large  measurement 
nonlinearities. 

The  nonlinearity  error  can  also  be  explained  by  graphically.  Figure  2.5  shows 
a  spherical  waveform  at  the  receiver  location.  The  uncertainty  orthogonal  to  the 
line-of-sight  from  transmitter  to  receiver  will  only  increase  the  range.  This  results 
in  the  range  to  be  under-estimated.  As  the  waveforms  become  more  spherical,  this 
error  becomes  more  substantial. 


Transmitter 


9 


line  of  sight 


Figure  2.5  Nonlinear  Elongation  of  Range 
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This  problem  is  not  unique  to  GPS  or  pseudolite  operation.  Widnall  [?]  referred 
to  this  problem  as  a  nonlinear  elongation  of  measured  range.  He  suggested  apply¬ 
ing  nonlinear  filtering  techniques  to  enlarge  the  region  of  convergence  of  a  Kalman 
filter.  He  stated  that,  if  the  nonlinearity  is  comparable  to  the  measurement  error, 
divergence  can  occur. 

2. 7  Summary 

This  chapter  has  provided  a  basic  overview  of  Kalman  filter  theory  including 
the  extended  Kalman  filters.  GPS  techniques  including  carrier-phase  differential 
algorithms  were  presented.  The  section  on  pseudolites  described  the  challenges  and 
issues  of  pseudolite  navigation. 
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III.  Methodology  and  Algorithm  Development 


3. 1  Overview 

This  chapter  begins  with  the  laying  out  the  simulation  design  and  setup.  It 
then  gives  a  description  of  the  truth  model  and  continues  with  the  generation  of  the 
measurement  and  measurement  model  errors.  Next,  the  floating-point  DGPS  filter  is 
described,  along  with  the  features  of  this  filter.  This  chapter  ends  with  a  description 
of  the  specific  carrier-phase  ambiguity  techniques  used  in  this  research. 

3.2  Overall  Simulation  Design 

This  research  involved  development  of  a  Kalman  filter-based  processing  algo¬ 
rithm  for  calculating  position  and  ambiguities  for  the  pseudolite  system,  as  well  as 
simulating  a  truth  model.  The  truth  model  was  not  a  traditional  truth  model,  in 
the  traditional  Kalman  filter  performance  evaluation  sense,  in  which  a  low  order 
filter  is  compared  to  a  higher  order  truth  model.  Rather,  the  truth  model  represents 
the  true  positions  of  the  transmitters  and  receivers  along  with  the  true  ambiguities, 
and  simulated  errors  in  the  simulated  measurements.  The  purpose  of  the  filter  algo¬ 
rithm  was  to  estimate  the  position  of  the  mobile  receiver  and  the  ambiguities  in  the 
carrier-phase  measurements. 

3.3  Truth  Model 

The  first  step  in  the  generation  of  the  truth  model  consisted  of  selecting  the 
coordinates  of  the  stationary  pseudolites  and  reference  receiver.  The  projected  coor¬ 
dinates  of  the  Inverted  GPS  Range  (IGR)  were  used  for  this  research  [16].  The  IGR 
is  a  GPS  modernization  field  test  program  which  is  currently  unfunded  at  the  time 
of  this  writing.  The  coordinates  reflect  realistic  pseudolite  locations  near  Holloman 
AFB,  NM.  Table  3.1  lists  the  positions  of  the  pseudolites  (numbered  1-10)  and  ref- 
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erence  receiver  while  their  positions  are  shown  in  Figure  3.1  relative  to  the  reference 


receiver. 


Table  3.1  Psenclolite  and  Reference  Receiver  Truth  Locations 


PRN  # 

Latitude  (deg) 

Longtitude  (deg) 

Ellipsoid  height  (nr) 

1 

33.50321 

-106.56055 

1433 

2 

33.54685 

-106.43650 

1630 

3 

33.67726 

-106.67454 

1447 

4 

33.39548 

-106.67449 

1541 

5 

33.66746 

-106.43772 

1602 

6 

33.78633 

-106.48492 

1562 

7 

33.82624 

-106.66540 

1574 

8 

33.64742 

-106.56846 

1426 

9 

33.72906 

-106.41564 

1656 

10 

33.72990 

-106.63744 

1442 

Ref  Rcvr 

33.65695 

-106.53696 

1424 

* 


-20  -15  -10 


O  Reference  RCVR 
*  Pseudolites 


-5  0  5 

Easting  (meters) 


10  15  20  25 


Figure  3.1  Psendolite  and  Reference  Receiver  Relative  Positions 
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The  second  step  in  the  truth  model  generation  was  the  selection  of  a  realistic 
flight  trajectory.  The  flight  trajectory  used  for  this  research  came  from  actual  flight 
test  data  collected  from  a  C-21  operated  by  the  746th  Test  Squadron  at  Holloman 
AFB.  The  simulation  is  constructed  so  only  the  positions,  i.e.  no  velocities  or  acceler¬ 
ations,  of  the  mobile  receiver  are  required  to  run  additional  flight  scenarios.  CHAPS 
is  normally  flown  on  a  C-21  from  Holloman  AFB,  NM,  so  the  flight  trajectory  used 
in  this  research  is  realistic  of  a  typical  reference  system  flight.  The  three-dimensional 
view  of  the  trajectory  with  projections  unto  each  axis  is  plotted  in  Figure  3.2. 


Figure  3.2  3-Dimensional  View  of  Trajectory  With  Projections 


Although  the  flight  trajectory  and  the  coordinates  of  the  pseudolites  and  refer¬ 
ence  receiver  are  located  near  Holloman  AFB,  New  Mexico,  they  are  not  co-located. 
The  coordinates  of  the  pseudolites  and  reference  receiver  were  adjusted  to  within 
range  of  the  mobile  receiver  flight  trajectory.  This  was  accomplished  by  first  con¬ 
verting  them  from  the  LLH  frame  to  an  ECEF  frame,  shifting  them  in  a  local  level 
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frame,  and  finally  converting  back  to  ECEF.  This  arrangement  was  necessary  to 
maintain  the  exact  relative  arrangement  of  the  psendolite  network. 

3-4  Measurement  and  Measurement  Model  Error  Generation 

Section  2.5.6  developed  the  donble-differenced  code  and  carrier-phase  measure¬ 
ment  equations  for  GPS  applications.  Those  equations  require  adaptation  for  use 
with  pseudolite  measurements.  The  ionospheric  error  terms  can  be  removed  because 
pseudolite  signals  will  not  travel  through  the  ionosphere.  The  resulting  equations 
for  pseudolite  applications  are 

AVpX,  =  AVrJj,  +  A  VT%  +  A  Vm%B  +  AVv»B  (3.1) 

for  the  code  and 

AV^b  =  A-‘(AVys  +  A  ,VI%  +  A  Vm%B  +  AV,^)  +  A  (3.2) 

for  the  carrier-phase.  GPS  navigation  is  affected  by  errors  in  the  predicted  motion  of 
the  satellites,  commonly  referred  to  as  ephemeris  or  orbital  errors.  This  error  is  not 
present  in  the  signal  itself,  but  occurs  when  the  receiver  uses  the  imprecise  satellite 
locations  for  range  calculations.  Pseudolites  have  a  corresponding  error  that  is  due 
to  the  imprecise  estimates  of  the  reference  receiver  and  pseudolite  locations.  These 
pseudolite  and  reference  receiver  position  errors,  along  with  the  measurement  noise, 
multipath,  and  tropospheric  delay  errors  terms  were  simulated  in  this  research. 

The  following  sections  describe  the  process  of  generating  the  measurement  and 
measurement  model  errors  used  in  this  research.  First  the  pseudolite  position  errors 
are  described,  followed  by  the  descriptions  of  measurement  noise,  multipath,  and 
tropospheric  delay.  Measurement  noise,  multipath,  and  tropospheric  delay  error 
were  added  to  the  true  ranges.  The  pseudolite  and  reference  receiver  position  errors 
were  added  to  the  true  positions,  and  those  positions  were  used  in  the  filter.  The 
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figures  in  each  section  show  plots  for  code  and  carrier-phase  errors  for  ten  pseudolites 
for  every  epoch  in  the  simulation.  These  plots  show  the  relative  magnitudes  and  time 
correlations  for  each  error  source.  After  the  errors  were  added  to  the  true  ranges, 
but  before  a  data  file  was  created  with  the  code  and  carrier-phase  measurements,  a 
maximum  range  limitation  was  implemented.  The  maximum  range  feature  simulated 
a  pseudolite  that  was  out  of  signal  range  from  the  mobile  receiver.  Pseudolite  prn  4 
was  31,698  meters  from  the  reference  receiver,  so  the  maximum  range  was  typically 
set  higher  than  that  at  32,000  meters.  This  resulted  in  7-10  satellites  visible  for  the 
mobile  receiver. 

3-4-1  Pseudolite  and  Reference  Receiver  Position  Errors.  The  imprecise 
estimated  positions  of  the  pseudolites  and  reference  receiver  affect  the  code  and 
carrier-phase  measurement  by  the  same  magnitude.  The  location  errors  of  the  pseu¬ 
dolites  and  reference  receiver  were  modeled  as  biases  with  a  zero-mean  Gaussian 
distribution.  Errors  were  created  in  an  East-North-Up  reference  frame  with  inde¬ 
pendent  horizontal  and  vertical  components  which  were  added  to  the  true  positions. 
The  horizontal  standard  deviation  was  set  to  1  cm  and  the  vertical  was  set  to  2  cm 
to  represent  the  expected  accuracies  of  precision  surveying. 

The  errors  due  to  inaccurate  positions  of  the  pseudolites  and  reference  receiver 
were  not  added  to  the  true  range,  but  instead  used  by  the  filer  in  the  measurement 
prediction  calculation.  Figure  3.3  shows  the  equivalent  effect  of  the  position  errors 
on  the  ranges  to  the  mobile  and  reference  receiver.  The  ranges  calculated  at  the 
mobile  receiver  include  only  the  effect  of  position  errors  at  the  pseudolites  while  the 
ranges  calculated  at  the  reference  receiver  include  both  the  position  errors  of  the 
pseudolites  and  the  errors  of  the  reference  receiver  itself. 

The  ten  lines  shown  on  Figure  3.3,  one  for  each  visible  pseudolite,  denote  the 
range  error.  The  lines  for  the  reference  receiver  range  errors  are  constants  because 
the  positions  errors  are  constant  during  a  simulation  run. 
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Figure  3.3  Position  Error  Effect  on  Ranges 

3-4-2  Measurement  Noise.  The  measurement  noise  was  modeled  as  zero- 
mean  white  Gaussian  noise.  Measurement  noise  is  considerably  smaller  for  carrier- 
phase  measurements  than  for  code  measurements,  and  it  was  modeled  with  a  55  cm 
standard  deviation  for  the  code  versus  a  3.5  mm  standard  deviation  for  the  carrier- 
phase  [28].  The  measurement  noise  was  modeled  as  uncorrelated  between  the  mobile 
and  reference  receiver. 

3-4-3  Multipath.  The  multipath  that  exists  in  a  pseudolite  system  is  caused 
by  the  environment  at  the  location  of  the  transmitters  and  receivers  [8,  12].  The 
reference  receiver  can  experience  more  time-invariant  multipath  than  the  mobile 
receiver  because  there  isn’t  any  relative  movement  between  the  reference  receiver 
and  the  pseudolites  .  These  pseudolite  phenomena  can  be  reduced  with  a  Multipath- 
Limiting- Antenna  on  both  the  transmitters  and  receivers  [8] .  For  the  purposes  of  this 
research,  the  multipath  for  the  received  signal  at  the  reference  receiver  was  modeled 
with  longer  time  correlations  than  the  multipath  for  the  mobile  receiver. 
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The  multipath  was  generated  as  a  first  order  Gauss  Markov  process  driven 
with  the  parameters  o  and  r,  which  are  themselves  modeled  as  first  order  gauss 
Markov  processes  (shown  below  as  FOGM(cr,  r).  A  minimum  threshold  was  set  for 
the  <7  and  r  values.  The  code  multipath  was  generated  as  1.5  times  the  summation 
of  source  1  and  source  2,  to  accurately  model  multipath  error.  These  equations  can 
be  summarized  by 


multipath  =  FOGM(a,  r)  (3.3) 

o=o0  +  FOGM(ctct,  Ta)  (3.4) 

r=To  +  FOGM(tt,  tt)  (3.5) 

(T  ^  <J min  (3-6) 

T  >  Tmin  (3.7) 


The  values  used  for  the  mobile  receiver  multipath  are  listed  in  Table  3.2  are  taken 
from. 

Table  3.2  Multipath  Parameter  Values 


Parameter 

Mob  Code 

(Source  1) 

Mob  Code 

(Source  2) 

Mobile 

Phase 

Ref  Code 

(Source  1) 

Ref  Code 

(Source  2) 

Reference 

Phase 

0o  (cm) 

10.0 

20.0 

0.19 

10.0 

20.0 

0.19 

r0  (sec) 

500 

25 

1000 

500 

25 

1000 

0V  (cm) 

4.0 

0.1 

0.038 

4.0 

0.1 

0.038 

ra  (sec) 

2000 

2000 

1500 

2000 

2000 

1500 

oT  (sec) 

200 

2 

400 

200 

2 

400 

rT  (sec) 

2000 

2000 

2000 

2000 

2000 

2000 

Omin  (cm) 

5.0 

1.0 

0.019 

5.0 

1.0 

0.019 

Tmin  (sec) 

100 

1 

0.019 

100 

1 

0.019 
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The  parameters  for  the  reference  receiver  multipath  were  the  same  for  the 
mobile  receiver  except  that  the  time  constants  were  tripled  to  simulate  stronger 
time  correlations.  Figure  3.4  depicts  a  typical  run  for  mobile  receiver  multipath 
error,  while  Figure  3.5  depicts  the  reference  receiver  multipath.  Both  plots  show  the 
code  multipath  on  the  top  subplot  and  the  carrier-phase  multipath  on  the  bottom 
subplot. 
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Figure  3.4 


Typical  Mobile  Receiver  Multipath  Error  for  All  Ten  Pseudolites 
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Figure  3.5  Typical  Reference  Receiver  Multipath  Error  for  All  Ten  Pseudolites 

By  comparing  Figures  3.5  and  3.4  it  is  clear  that  the  multipath  simulated  at 
the  reference  receiver  had  strong  time  correlations  between  epochs.  Recall  that  this 
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was  the  purpose  of  tripling  the  time  constants  for  the  reference  receiver  tropospheric 
truth  model. 


3-4-4  Tropospheric  Delay.  The  truth  model  for  generating  the  tropospheric 
delay  was  taken  from  reference  [3,  33]  which  was  a  function  of  temperature,  atmo¬ 
spheric  pressure,  relative  humidity,  elevation  angle,  and  range.  The  atmospheric 
parameters  are  the  ones  taken  at  the  reference  receiver.  The  tropospheric  delay 
calculation  for  the  mobile  receiver  is  defined  as 


taplARu,  A/lu) 

77.6 Ps  x  (42700  -  hs)  x  10“6 
5TsAhu 

Ns  x  (13000  -  h8)  x  10“6 
5  A  hu 


^7~v,dry  +  A TVjWet  (A Tv, dry  T  ATv,wet) 


1  - 


sin  (elu) 

AHapl 
42700  -  ht 

AHapl 
13000  -  h. 


Ah ,, 


Rv 


-  1- 


A  Hapl  +  A  hv 
42700  -  hs 

A  Kapl  +  Ahv 
13000  -  hq 


Rv 


Rh 


(3.8) 


where  the  variables  are  defined  as 


TAPL,u 

R-u 

Ahu 


At , 


v,dry 


At. , 


v,wet 


f>] 


A  Hapl 


hs 


Ts 

Ns 


=  tropospheric  delay  for  mobile  receiver  (meters) 

=  slant  range  between  the  pseudolite  and  user  (meters) 

=  the  height  of  the  user  above  the  pseudolite  (meters) 

=  differential  vertical  dry  delay  (meters) 

=  differential  vertical  wet  delay  (meters) 

=  elevation  angle  in  radians 

=  difference  in  height  between  pseudolites  and  reference  receiver 
=  height  of  reference  receiver 
=  surface  pressure  (millibars) 

=  surface  temperature  (Kelvins) 

=  surface  refractivity 
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Equation  3.8  is  adapted  for  reference  receiver  tropospheric  delay  calculation 


and  given  by 


tapl,r{Rr,  AIir)  — 


+ 


77.6 Ps  x  (42700  -  ha)  x  10~6 
5TsA1iapl 

Ns  x  (13000  -  h8)  x  10~6 


1  -  1  - 


Ah 


APL 


42700  -  h. 


5  A  Iiapl 


1  -  1  - 


Ah 


APL 


13000  -  h. 


Rr 


Rr  (3.9) 


Both  Equations  3.8  and  3.9  are  valid  for  positive  and  negative  elevation  angles  but 
indeterminate  for  zero  elevation  angles  [3,  33].  The  reference  did  develop  equations 
for  zero  elevation  angles,  but  they  were  not  implemented  in  this  research,  because 
the  placement  of  the  reference  receivers  and  mobile  receiver  trajectory  did  not  result 
in  zero  elevation  angles. 


3.5  Floating-Point  Differential  Pseudolite  Kalman  Filter 

The  floating-point  differential  pseudolite  Kalman  filter  used  in  this  research 
is  a  post-processed  algorithm.  It  is  a  modified  version  of  the  filter  developed  in 
reference  [27],  which  is  adapted  for  pseudolite  navigation.  The  double-difference 
operation  is  applied  to  both  the  code  and  carrier-phase  measurements,  allowing 
the  removal  of  the  states  modelling  clock  error.  This  section  presents  the  baseline 
filter  development,  and  the  modifications  are  presented  in  the  next  section.  The 
baseline  filter  calculates  position,  velocity,  acceleration,  and  carrier-phase  ambiguity 
estimates  for  the  mobile  receiver.  The  objective  of  the  filter  is  to  produce  carrier- 
phase  ambiguity  estimates  and  associated  covariances  that  will  be  processed  through 
ambiguity  resolution  techniques  to  produce  the  fixed-integer  results. 

Before  the  filter  is  run  on  the  data,  a  pre-processing  step  is  conducted  to 
determine  the  number  of  visible  pseudolites,  a  vector  of  visible  pseudolite  prns,  the 
base  pseudolite  for  double-difference  operation,  and  the  non-base  pseudolite  prns  for 
each  epoch  data. 


3-10 


3.5. 1  Differential  Pseudolite  Model  Equations.  A  First  Order  Gauss  Markov 
Acceleration  (FOGMA)  model  was  used  to  define  the  three  positions,  three  velocities, 
and  three  accelerations  states  of  the  floating-point  differential  pseudolite  Kalman  fil¬ 
ter.  The  remaining  states  consisted  of  (n  —  1)  carrier-phase  ambiguity  states,  where 
n  is  the  number  of  pseudolites  in  view  at  that  epoch. 

The  positions,  given  in  the  ECEF  coordinate  frame,  are  modeled  as  the  time 
derivatives  of  the  velocities,  and  the  velocities  are  modeled  as  the  time  derivatives 
of  the  accelerations  by: 


X\  =  X4  X4  —  X7 

x-2  =  x5  x5  =  x8 

X-3  =  X6  X6  =  Xg 


(3.10) 


The  position  and  velocity  states  are  completely  determined  by  other  states,  and 
they  do  no  include  any  direct  driving  noise.  The  acceleration  states  are  modeled  as 
first-order  Gauss- Markov  processes  by 

X7  =  (-1  /Ta)x7  +  wai(t) 

x8  =  (-1  /Ta)x8  +  wa2(t)  (3.11) 

Xg  =  (~1/Ta)xg  +  Wa3(t) 

with  associated  dynamic  driving  noise  processes  given  by 

E  {wai  ( t)wai  (t  +  t)}  —  E  {wa2  ( t)wa2  (t  +  r)}  =  E  {wa3 (t) waz  (t  +  r) } 

=  ^M=?a. 5(r)  (3.12) 

a 

The  correlation  time,  Ta,  and  variance  (or  mean  square  value),  af  for  the  accel¬ 
erations  are  determined  based  on  the  anticipated  acceleration  maneuvers  and  time 
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correlations.  The  value  Ta  was  set  to  3  seconds  to  account  for  relatively  short  accel¬ 
eration  maneuvers  in  a  C-12,  and  cra  was  set  to  15  m/s 2  to  handle  the  ’’worst  case” 
acceleration.  With  these  values  qa  is  calculated  to  be  150  m2  /  sec5. 

The  floating-point  filter  can  accept  either  a  single  or  widelane  frequency.  In 
either  case,  the  double  difference  is  performed  with  the  ambiguity  terms  still  main¬ 
taining  an  integer  nature.  The  ambiguities  were  modeled  as  random  walks  rather 
than  constant  biases  to  ensure  that,  if  the  filter  converged  to  an  incorrect  value, 
it  could  correct  itself  (i.e.,  the  gain  of  the  .  The  cycle  ambiguities  were  modeled 
in  an  additional  (n  —  1)  states  after  the  initial  9  states  for  position,  velocity,  and 
acceleration.  The  double-differenced  carrier-phase  ambiguities  are  defined  by: 


xio  —  WvAAU-2 


Xu  —  Wv  AN1-3 


X(8 +n)  ~  WVAJV1-"  (3.13) 

where  PRN  1  is  given  as  the  base  and  n  represents  the  total  number  of  pseudolites 
visible. 

The  process  noise  is  given  as 

E{wVANbi(t)wVANu(t  +  r)}  =  qN8(T) 
qN  =  1.1  x  10-5  (cycles2/ sec) 

The  value  of  qA  will  yield  an  increase  of  approximately  0.2  cycles  in  the  ambiguity 
standard  deviation  over  a  1  hour  period  [27].  This  will  allow  the  filter  to  correct 
itself  if  it  converged  to  the  incorrect  value. 
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The  state  vector  for  the  floating-point  Kalman  filter  is  defined  as 

x  =  [XYZXYZXYZ  VAN1'2  . . .  VAN1'11]  T  (3.14) 

where 

X\  =  X  =  ECEF  X  position  (m) 
x2  —  Y  —  ECEF  Y  position  (m) 
x3  =  Z  =  ECEF  Z  position  (m) 

X4  =  X  =  ECEF  X  velocity  (m/s) 

X5  —  Y  —  ECEF  Y  velocity  (m/s) 
xe  =  Z  =  ECEF  Z  velocity  (m/s) 

£7  =  X  =  ECEF  X  acceleration  (m/s2) 
x§  —  Y  —  ECEF  Y  acceleration  (m/s2) 
xg  —  Z  —  ECEF  Z  acceleration  (m/s2) 

£10  =  VAN1'2  =  double-differenced  phase  ambiguity  (cycles) 

£11  =  VAN1'3  =  double-differenced  phase  ambiguity  (cycles) 

£g_i_(n_i)  =  VAN1'”  =  double-differenced  phase  ambiguity  (cycles) 

The  differential  equation  is  similar  to  Equation  2.1  and  is  represented  by 

x(t)  =  F  (t)x(t)  +  G  (t)w(t)  (3.15) 
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which  expands  to: 


Xi 

0 

0 

0 

1 

0 

0 

0 
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0 

0  ••• 

0 

X\ 

0 

X2 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0  ••• 
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X2 

0 

x3 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0  ••• 

0 

x3 

0 

X  4 

0 

0 

0 

0 
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0 

1 

0 

0 

0 

0  ••• 

0 

X4 

0 

X5 

0 

0 

0 
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0 

0 

1 

0 

0 

0  ••• 

0 

x5 

0 

Xq 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0  ••• 

0 

x6 

0 

X7 

= 

0 

0 

0 

0 

0 

0 

-1  /Ta 

0 

0 

0 

0 

0 

x7 

+ 

wa  1 

Xs 

0 

0 

0 

0 

0 

0 

0 

—  1  /Ta 

0 

0 

0  ••• 

0 

x8 

Wa2 

Xg 

0 

0 

0 

0 

0 

0 

0 

0 

-1  /Ta 

0 

0  ••• 
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Xg 

Wa3 

Xio 
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Xio 

Nbi 

in 
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0 
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0 
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0 

0  ••• 

0 

Xu 

w\7ANbi 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0  0 

0 

Nbi 

(3.16) 


In  Equation  3.16,  the  variable  Ta  represents  the  FOGMA  acceleration  time  constant. 
The  G  (t)  matrix  is  defined  to  be  an  identity  matrix  for  this  research.  The  dynamics 
driving  noise  Q  is  defined  by: 

E{w(t)w7  (t  +  r)}  =  Q(t)5(r)  (3-17) 
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The  matrix  Q  is  represented  by: 
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0  0 
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The  acceleration  mean  squared  value,  time  constant,  and  acceleration  noise,  and  the 
phase  ambiguity  noise  values,  which  were  previously  justified,  are  summarized  in 
Table  3.3. 

Table  3.3  Floating-Point  Filter  Dynamics  Driving  Noise  Values 


Term 

Definition 

Value 

Mean  squared  value 

(12.25  m/sec2)2 

T 

Acceleration  time  constant 

3  seconds 

Qa 

Acceleration  noise 

100  m2  j sec 5 

Qn 

Phase  ambiguity  noise 

1.1  x  10“4  cycles2 / sec 

The  initial  conditions  for  the  position  states  were  set  to  the  true  value  at  the 
first  epoch  with  an  additive  zero-mean  error  term  that  had  a  Gaussian  distribution 
and  a  standard  deviation  of  5  meters.  The  velocity  and  acceleration  initializations 
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were  set  to  zero  with  the  covariance 

initialization  matrix  defined  by: 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

al 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

a'2y 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

P(«o)  = 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

2 

°VA  Nbi 
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ffVAJVH  •  •  • 
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0  0 

2 

°VA  Nbi 

(3 

The  initial  covariance  values  used  in  this  research  are  given  in  Table  3.4  where  the 
number  of  visible  pseudolites  in  the  first  epoch  is  given  as  n. 

Table  3.4  Floating-Point  Filter  Initial  Covariance  Values 


Term 

Definition 

Value 

2 

® x,y,z 

Position  state  variance 

(100  mf 

2 

ax,y,z 

Velocity  State  variance 

O 

o 

to 

2 

® x,y,z 

Acceleration  state  variance 

(20  m/s2)2 

2 

aVANbi 

Phase  ambiguity  variance 

(y  cycles )2 

3.5.2  Differential  Pseudolite  Measurement  Model.  The  floating-point  dif¬ 
ferential  pseudolite  Kalman  filter  uses  a  nonlinear  measurement  model  which  consists 
of  double-differenced  code  and  phase  measurements  resulting  in  a  2 (n  —  1)  measure- 
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ment  vector  (where  n  is  the  number  of  pseudolites  in  view)  defined  by 

z(ti)  =  [VA p1'2  ■  ■  ■  VAp1~nVA(/)1~2  ■  ■  ■  \7A(j)1-n}T  (3.20) 

From  Chapter  2,  Equation  2.23  described  the  nonlinear  measurement  model  for  an 
extended  Kalman  filter  which  is  in  the  form: 

z  (U)  =  h[x(C),  tj]  +  v(C)  (3.21) 

It  must  be  linearized  before  it  is  used  in  the  gain  and  covariance  calculations  of  the 
extended  Kalman  filter.  Recall  from  Chapter  2  that  thee  partial  derivative  matrix 
H  was  defined  as 

(3.22) 

x=x(t“ ) 

which  is  an  m  x  s  matrix,  where  m  is  the  number  of  measurements  and  s  is  the 
number  of  states.  Each  row  corresponds  to  a  single  measurement  and  is  defined  as 

<9h[x,  tj] 

<9x 

x=x(ti  ) 

Recall  from  Equation  2.57  that  the  double-differenced  carrier-phase  measurement  is 
given  by 

A  V<t>kJs  =  A'1  (A  Vr%  +  A  VT$JB  +  A  Vm%  +  AVv%)  +  A' V  NkJB  (3.24) 

The  A  term  in  Equation  3.24  is  the  carrier  wavelength  and  will  depend  on  whether 
a  single  or  widelane  frequency  is  used. 

When  the  double-differenced  range  term  is  expanded  and  the  measurement 
errors  are  combined  the  carrier-phase  equation  is  expressed  as 

AV0ab  =  j  [- r°B  -  r\  -  ( rkB  -  rkA)\  +  A WN3AB  +  vAVlf)  (3.25) 


<9h[x,  U] 

<9h[x,C] 

<9h[x,C] 

(9xi 

x=x(t . 

(9x2 

) 

x=x(t .  ) 

dxs 

x=x(i.  ) 

(3.23) 


[x,ti 

H[C;x(t,.  )]  = 


(9x 
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In  the  preceding  equation,  the  term  vav</>  is  modeled  as  a  white  noise,  and  it  repre¬ 
sents  the  combination  of  the  doubled-differenced  measurement  noise,  multipath,  and 
tropospheric  delay.  It  is  important  to  note  that  the  tropospheric  delay  terms  are  the 
residual  error  after  a  tropospheric  model  has  been  applied.  When  the  range  terms 
in  Equation  3.25  are  further  expanded  and  expressed  in  terms  of  state  variables,  the 
equation  becomes 

A V =  j  [(x3  -  x3)2  +  ( y 3  -  x2f  +  (zJ  -  x3)2] 1/2 

-  [( xk  -  Xl)2  +  (: yk  -  x2)2  +  (zk  -  x3)2} 1/2  (3.26) 

+  j  {rkA  -  r3A)  +  AVN3AkB  +  wAv<i» 

where  x3,k ,y3,k,  z3,k  represent  the  estimate  of  the  pseudolites  indexed  by  j  and  k. 
Recall  the  y  {rB  —  rJR }  term  represents  the  mobile  receiver  and  is  expanded,  but  the 
X  {rA  ~  ryi}  term  represents  the  reference  receiver  so  it  is  not  a  function  of  state 
variables,  and  thus  not  expanded. 

The  partial  derivatives  for  each  row  of  the  double-differenced  carrier-phase 
measurements  are  given  as 

d 7i[x,  U\  I 


^  l  [(xj-xi)2+(yi-x2)2+(zj-X3)2]1/2 


^(xk-xi)2-\-(yk—X2)2-\-(zk—xs)2^  ^ 

=  Hei~ei} 


(3.27) 


d  /i[x,  ti 


dx2 


X=x(ti  ) 


_  _ y:  —x  2 _ 

^  \  [(xi -xi)2+(yi —X2)2+(zi —X3)2]132 

.1  / _ yk-x  3 _ 

^  \  \L(xk—x\)2+(yk—X2)'2+(zk—xz) 2]1^“ 

=  lie  2-4} 


(3.28) 
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1 

A 


_ z3  -3:3 _ 

[{x3  —Xi)2+(y3  —X2)2+(z3  -X3)2]1/2 


dh  [x,  ti 


dx3 


x=x(ti  ) 


! 

A 


[(xfc  —  Xl)2  +  (j/fe—  X2)2+(zk- CC3)2] 

=  I  -  4} 


1/2 


(3.29) 


where 


d  /i[x,  h; 


(9x 


bi 


=  1 


x=x(ti  ) 


(3.30) 


fv?  = 

mob 


i  1 
e\  eJ2 


(3.31) 


is  the  unit  line-of-sight  vectors  pointing  from  the  mobile  receiver  to  psendolite  j. 

When  these  individual  partial  derivatives  are  combined  they  represent  one  row 
of  the  H  matrix  as 


Wk 


A 


0  0  0  0  0  0 


(3.32) 


where  j(eJmob  —  ekmob)  represents  the  scaled  difference  vector  between  two  unit  line- 
of-site  vectors  from  the  mobile  receiver  to  pseudolite  "  j"  and  the  mobile  receiver  to 
pseudolite  "k" .  The  ”1”  is  placed  in  the  column  for  the  appropriate  ambiguity  state. 

The  corresponding  rows  for  the  double  differenced  code  measurements  are  the 
same  values  after  the  j  term  is  removed  and  the  ”  1”  for  the  ambiguity  state  values 
is  dropped: 

II  ‘  =  [(<ot  -  elt)  0  0  0  0  0  0  ...  0]  (3.33) 
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The  entire  measurement  matrix  H  is  an  9+(n-l)  by  2(n-l)  matrix  defined  by 


(e1  -  e2) 

0 

0 

0 

0 

0 

0 

0 

0  •• 

— 

o 

(e1  -  e3) 

0 

0 

0 

0 

0 

0 

0 

0 

•  0 

(e1  —  en) 

0 

0 

0 

0 

0 

0 

0 

0  •• 

■  0 

K^-e2) 

0 

0 

0 

0 

0 

0 

1 

0 

■  0 

!(e\-e3) 

0 

0 

0 

0 

0 

0 

0 

1  •• 

•  0 

.  0 

-f  (e1  —  en) 

0 

0 

0 

0 

0 

0 

0 

0  •• 

•  1 

(3.34) 


where  b  is  the  base  pseudolite  for  n  pseudolites  where  i  =  1  •  •  •  n  ,  with  i  ^  base 
pseudolite. 

The  measurement  covariance  matrix  R  is  defined  by 


£'{v(^)vt(^)} 


R iti)  for  ti  ^  tj 
0  for  ti  ^  tj 


(3.35) 


is  required  by  the  filter.  This  matrix  can  be  broken  up  into  4  different  types  of 
covariance  terms. 


•  Case  1:  Variance  of  code  measurement  errors 

•  Case  2:  Variance  of  phase  measurement  errors 

•  Case  3:  Covariance  of  two  different  code  measurement  errors 

•  Case  4:  Covariance  of  two  different  phase  measurement  errors 
The  full  R  matrix  can  be  partitioned  into  four  sections  represented  by 

■Encode  0 

0  R. phase 
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where  the  code  variance  and  covariances  denoted  by  cases  1  and  3,  respectively, 
are  placed  in  the  upper  left  corner.  The  phase  variances  and  covariances  denoted 
by  cases  2  and  4  are  located  in  the  lower  right  corner.  The  upper  right  and  lower 
left  corners  represent  the  cross-covariance  of  a  code  and  phase  measurement.  These 
values  were  assumed  to  be  sufficiently  small  to  ignore  them. 

The  variances  for  the  code  measurement  are  a  combination  of  residual  tro¬ 
pospheric  error  and  the  non-tropospheric  components  (transmitter  location  error, 
multipath,  and  measurement  noise).  For  this  research,  the  non-tropospheric  er¬ 
ror  was  assumed  to  be  uncorrelated  between  measurements,  which  means  that  the 
double-differenced  standard  deviation  is  a  factor  of  two  greater  than  the  observation 
standard  deviation.  The  double-differenced  standard  deviations  for  the  tropospheric 
and  non-tropospheric  errors  were  3.2  meters  and  0.07  meters,  respectively,  which 
resulted  in  a  total  standard  deviation  of  3.2008  meters.  It  is  important  to  note  that 
the  tropospheric  contribution  to  the  total  standard  deviation  is  sufficiently  small 
to  ignore,  but  is  included  for  completeness.  The  covariances  of  two  different  code 
measurements  were  set  as  one  half  of  the  code  variances,  because  half  of  the  mea¬ 
surements  are  in  common  due  to  double  differencing.  The  following  matrix  defines 
the  code  segment  of  the  R  matrix. 


R, 


■code 


^VAp4-?  ,VAp*-?  ^VAp*-?  ,VAp4fc 

^VAp^  ,VApifc  ^VAp^  ,VAp^ 

rVApii,VAp!fe 


^VApO',VApifc 

1  VApVjVAplfc 
1  VAp'J,VAp8fc  ^  VAp'J, VAp'J 


(3.36) 


The  phase  variances  and  covariances  were  developed  in  a  similar  manner  to  the 
code  values.  The  double-differenced  tropospheric  and  non-tropospheric  standard  de¬ 
viations  are  0.0812  cycles  and  0.0464  cycles,  respectively,  with  the  resulting  total 
standard  deviation  of  0.0935  cycles.  The  phase  variance  was  calculated  to  be  0.0087 
square  cycles  and  the  covariances  0.00435  square  cycles.  The  following  matrix  dis- 
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plays  the  phase  component  of  the  full  R  matrix. 


R 


'phase 


^VA^’,VA^'  ^VA^’,VA^fc 
^VA^’,VA(/)ifc  ^VA^,VA  </>*•? 


^VA^,VA^fc 

^VA^,VA 

^VA^,VA(/)*fc  ^VA0^',VA  (j>^ 


(3.37) 


The  values  for  the  the  R  matrix  are  shown  in  Table  3.5. 


Table  3.5  Measurement  Covariance  Values 


Term 

Definition 

Value 

^VA  pli ,  V  A  p*-? 

Double-differenced  code  variance  error 

10.24  m 2 

rvA/yqvApife 

Double-differenced  code  covariance  error 

5.12  m2 

r  VA0O,VA</>O 

Double-differenced  carrier-phase  variance  error 

.0087  cycles 2 

^VA^,VA^fc 

Double-differenced  carrier-phase  covariance  error 

.00435  cycles 2 

3.5.3  Discrete-Time  Models.  The  linear  stochastic  differential  equations 
must  be  converted  to  be  implemented  on  a  digital  computer.  This  requires  the 
formulation  of  the  linear  stochastic  difference  equations  to  describe  the  equivalent 
discrete-time  system  model  [19],  which  is  in  the  form 


x(4+i)  =  3*(4+i,4)x(4)  +wd 


(3.38) 


where 


E{wd}  =  0 

4{wd(4)wJ(4)}  =  Qd  (3.39) 

E{™d(tj)wd(tk)}  =  0,  tj  ^  tk 
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Equations  3.38  through  3.40  were  derived  in  Chapter  2.  The  discrete-time  state 
transition  matrix  G>{tk+i,tk)  is  defined  as: 


$(4+i,4)  =  $(At)  =  eFA*  (3.40) 

where  At  =  tk+ 1  —  tk  which  results  in  the  matrix 

100  At  0  0  AO  0  0  0  ■  •  ■  0 

0100  AtO  0  AO  0  0  ■  ■  ■  0 
0010  0  At  0  0  A  0  0  ■  ■  ■  0 

0001  0  0  50000  •••0 

0000  1  0  05000  •••0 

0000  0  1  00500  •••0 

$(4+i,4)=  0000  0  0  CO  0  0  0  ■  ■  ■  0  (3.41) 

0000  0  0  0  C  0  0  0  •••  0 

0000  0  0  0  0  C  0  0  0 

0000  0  0  00010  •••0 

0000  0  0  00001  •••0 

;  ;  ;  ;  ;  ;  o 

0000  0  0  000000  1 

where 


A  =  T*(e~At/Ta  -  1)  +  TaAt 
B  =  Ta{l  -  e~At/Ta) 

C  =  (e”A</T“) 
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The  discrete  dynamics  driving  noise  is  given  by 


rtk+i 

Q d(tk)  =  /  $(4+i  t)G(t)Q(t)Gt(t)$t(4+iiT)</t 

Jtk 

which  expands  to 

40040  0  G0  0  00  •••0 

0  D  0  0  E  0  0  GO  00  •••0 

0  0  4  0  0  4  0  0  G  0  0  •••  0 

E00K00L0  0  0  0  •  •  •  0 

0400  4  0  0  L  0  0  0  •••0 

0040  0  /1  0  0  L  0  0  ■  ■  ■  0 

Qd=  GO  0  LO  0  MO  0  0  0  •••0 

0G00L00  MO  00  ■■■0 

0  0  GO  0  LO  0  MO  0  ■■■0 

0000000  0  0  iVO-'-O 

0000000  0  0  0  N  •••  0 

000000000000  N 


(3.42) 


(3.43) 


3-24 


where  these  values  were  taken  from  reference  [27]  and  defined  as 


D  =  ir„59„(  1  -  e-2At/ra)  +  T‘rhM(1  _  2~MIT‘)  -  T3<ja(A«)2  +  ir„2?a(A«)3 
E  =  T*M  +  i)  +  r34«Af(e-A‘2T'  -  1)  +  ^.(Ai)2 

G  =  tr3?a(l  -  e-^‘/T-)  -  T*qal\te-A,/T‘ 

I<  =  +  2^-3) 

^  J-  a 

L  =  ~T^qa(-e~2At/Ta  +  2e~At/Ta  -  1) 

M  =  -Iraga(-e-2Ai/T“  -  1) 

N  =  qN^t 


3. 6  Floating-Point  Filter  Features 

The  floating-point  filter  included  a  pre-filtering  step,  real-data  considerations, 
a  tropospheric  model,  and  adaptations  to  improve  upon  the  performance  of  the 
baseline  filter.  These  adaptations  consisted  of  optimal  smoothing  techniques,  second 
order  filtering,  weighted  measurement  covariance  matrix,  and  estimating  errors  in 
the  tropospheric  model. 

3.6.1  Pre-filter.  A  pre-filter  function  was  implemented  to  determine  the 
number  of  available  PRNs,  vector  of  available  PRNs,  and  base  PRN  for  double 
difference  operations.  The  base  PRN  was  initially  chosen  from  a  vector  of  prns  that 
were  in  view  at  the  first  epoch.  The  PRN  that  stayed  in  view  the  longest  was  chosen. 
If  that  PRN  did  not  stay  in  view  the  entire  data  set,  the  process  was  repeated  in  a 
similiar  manner.  The  process  then  defined  a  vector  of  PRNs  that  were  in  view  at  the 
second-to-last  epoch  the  initial  base  PRN  went  out  of  view.  The  PRN  that  stayed  in 
view  the  longest  was  chosen  as  the  base,  and  if  it  did  not  stay  in  view  until  the  end 
of  the  data  set,  the  process  repeated  itself.  The  pre-filter  function  then  generates  a 
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pseudolite  visibility  plot  denoting  the  base  PRN  with  a  thick  line.  Figure  3.6  shows 
a  typical  pseudolite  visibility  plot. 


08:04  08:06  08:09  08:11  08:13  08:16  08:18 

GPS  Week  Seconds\Local  Time 


Figure  3.6  Pseudolite  Visibility  Plot 


3.6.2  Real  Data  Considerations.  Although  only  simulated  data  was  used  in 
this  research,  every  attempt  was  made  to  make  it  as  realistic  as  possible.  The  filter 
was  designed  to  handle  these  real  data  considerations.  The  algorithm  developed  in 
this  thesis  included  the  ability  to  handle  pseudolites  going  out  of  view,  pseudolites 
coming  into  view,  and  a  change  in  the  base  double  difference  PRN. 

In  a  real-world  system,  pseudolites  will  go  out  of  view.  When  this  happened, 
the  filter  eliminated  the  appropriate  state  estimate  and  the  rows  and  columns  asso¬ 
ciated  with  this  pseudolite  prn  in  the  covariance  matrix.  For  example,  if  pseudolites 
1  through  5  were  visible  with  prn  1  as  the  double  difference  base,  sample  ambiguity 
state  values  could  be 

xio  =  VAW12  =  -2837.24 
xn  =  WAN13  =  10314.35 
x12  =  WAN14  =  -563.10 
x13  =  VAN15  =  124.73 
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These  four  ambiguity  states  would  be  in  states  10  through  13,  because  the  first 
9  states  represent  the  3  position,  3  velocity,  and  3  acceleration  states.  If  pseudolite 
4  went  out  of  view,  the  new  ambiguity  vector  would  be 

x10  =  WAN12  =  -2837.24 
xn  =  VAN13  =  10314.35 
X12  =  VAN15  =  124.73 

The  covariance  matrix  P  also  requires  adjusting.  In  this  example,  it  would 
have  initially  been  a  13  by  13  matrix.  By  way  of  example,  consider  a  case  in  which 
the  4-by-4  partition  of  the  ambiguity  variances  and  covariances  is 

0.0063449  0.0037111  0.004901  0.00066232 

0.0037111  0.059193  0.0038798  0.0012557 

P  = 

0.004901  0.0038798  0.0064154  0.001028 

0.00066232  0.0012557  0.001028  0.0032174 

In  order  to  remove  the  VA1V14  state  the  second  to  last  row  and  second  to  last  column 
are  eliminated.  As  a  result,  the  covariance  becomes 

0.0063449  0.0037111  0.00066232 

P  =  0.0037111  0.059193  0.0012557 

0.00066232  0.0012557  0.0032174 

In  a  real-world  system,  pseudolites  will  also  come  into  view.  When  a  pseudolite  comes 
into  view,  both  the  length  of  the  state  vector  and  the  dimensions  of  the  covariance 
will  have  to  increase  by  one.  The  first  step  is  to  estimate  the  new  ambiguity  term  in 
a  similar  fashion  to  the  way  the  ambiguity  states  were  initialized  at  the  first  epoch; 
see  Section  3.4.1.  The  covariance  matrix  variance  term  is  set  to  the  same  ambiguity 
variance  term  used  in  the  first  epoch,  (y cycles)2 ,  while  the  off-diagonal  terms  for 
the  row  and  column  are  set  to  zero.  Taking  the  previous  example  with  pseudolite  4 
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coming  back  into  view,  including  the  estimate  of  -564.1  for  VAA14  determined  by 
the  code,  would  be 


XlO 

=  VA  A12  = 

-2837.24 

Xu 

=  VAN 13  = 

10314.35 

Xl2 

=  VAA14  = 

-564.1 

Xl3 

=  VAA15  = 

124.73 

with  the  covariance  given  by 

0.0063449  0.0037111 

0.0037111  0.059193 

P  = 

0  0 
0.00066232  0.0012557 


0  0.00066232 

0  0.0012557 

3350  0 

0  0.0032174 


The  value  of  3350  for  the  new  ambiguity  state  variance  assume  a  widelane  wave¬ 
length. 

The  base  pseudolite  PRN  used  for  double  difference  operations  cannot  be  as¬ 
sumed  constant  over  the  entire  data  set.  This  thesis  included  a  function  that  trans¬ 
lated  the  ambiguities  from  that  last  epoch  to  the  current.  A  transformation  matrix 
was  formed,  based  on  the  available  prns  and  base  prn  at  both  the  last  epoch  and 
current  epoch.  The  double  differenced  ambiguities  are  combinations  of  single  dif¬ 
ferenced  ambiguities  which  makes  this  transformation  possible.  No  information  is 
lost  in  the  translation:  the  new  state  and  covariance  are  just  different  measurement 
combinations.  This  function  must  be  performed  before  the  ambiguity  states  and 
covariance  matrices  are  adjusted  for  lost  or  gained  pseudolites. 

The  relation  is  given  as: 


^ new  Tx0/d 

(3.44) 

■p  _  T'P 

r  new  r  old 

(3.45) 
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Consider  using  the  previous  example,  when  all  five  pseudolites  were  in  view  and  the 
current  base  PRN  was  going  out  of  view,  the  translation  matrix  T  to  switch  the 
double  difference  base  PRN  from  1  to  3  would  be 


10000000000  00 
01000000000  00 
00100000000  00 
00010000000  00 
00001000000  00 
00000100000  00 
T  =  00000010000  00  (3.46) 

00000001000  00 
00000000100  00 
0000000000-100 
0000000001-100 
0000000000-110 
0000000000-101 

Note  that  the  top  left  corner  is  a  9  by  9  identity  matrix,  which  preserves  the  position, 
velocity,  and  acceleration  estimates.  The  bottom  righthand  corner  values  (which  are 
bolded)  represent  the  portion  of  T  that  rearranges  the  ambiguity  estimates.  The 
same  transformation  matrix.  T,  is  applied  to  the  covariance  according  to  Equation 
3.45.  When  T  is  multiplied  by  x0^  it  forms 


VA1V31 

v  new 

=  -VAiV« 

VA  N32 

v  new 

=  -VAJV“ 

+VAAPfd 

VA  A34 

new 

=  -VAA 

+vaa^ 

VA  A35 

v  new 

= -VAAS 

+VAA33 
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3.6.3  Tropospheric  Model.  The  tropospheric  model  used  in  the  filter  was 
the  same  as  that  used  to  generate  the  tropospheric  delay  in  the  truth  simulation. 
There  are  two  differences  between  the  truth  and  filter  model  for  the  tropospheric 
delay  calculation. 

The  first  source  of  residual  tropospheric  error  was  due  to  the  atmospheric  values 
used  in  each  model.  The  truth  model  used  true  values  for  atmospheric  pressure,  tem¬ 
perature,  and  relative  humidity.  Errors  were  added  to  the  true  atmospheric  values, 
which  were  then  used  in  the  filter  model.  These  errors  were  modeled  as  zero-mean 
Gaussian  random  biases  with  adjustable  standard  deviations.  The  standard  devia¬ 
tions  were  set  as  4  percent  for  relative  humidity,  1  degree  Kelvin  for  temperature, 
and  3  millibars  for  the  pressure. 

The  second  error  source  was  from  the  difference  in  using  the  estimated  versus 
the  true  positions  of  the  receivers  and  transmitters.  The  truth  model  used  the  true 
positions  while  the  filter  model  used  the  estimated  positions.  This  difference  resulted 
in  very  small,  essentially  insignificant,  errors  for  both  the  slant  ranges  and  elevation 
angles.  This  difference  was  very  small  in  magnitude  because  the  position  errors  were 
only  a  few  centimeters  while  the  ranges  were  measured  in  kilometers. 

3.6.4  Optimal  Smoothing.  The  optimal  smoothing  algorithm  presented  in 
Section  2.4  was  implemented  to  improve  the  accuracy  of  the  state  estimates  while 
decreasing  the  size  of  the  covariance  values.  The  algorithm  required  modifications 
when  pseudolites  were  allowed  to  come  in  and  out  of  view  and  the  base  PRN  for 
double  difference  operations  changed  between  epochs.  The  equations  give  in  Chapter 
2  are  restated  here  as 
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*(*/!</)  =*(</) 

=  *(i+)  +  A(fj)[x(ti+1|f/)  -  x(i7+1)] 

A(U)  =P(«,+)#T(i,+i,t,)P-'(<-+i) 

P(t,M  =P(tp 

P(t,M  =  P(«,+)  +  A((,)[P((,+1|(/)  -  P((-+1)]AT((,) 


which  are  run  backward  from  the  final  time  tf  after  a  first  forward  pass  is  made 
through  the  data  with  the  extended  Kalman  filter. 

These  fixed-interval  smoother  equations  obtained  by  Meditch  [19]  will  fail  when 
the  number  of  states  from  one  epoch  to  another  is  not  consistent,  or  the  quantities 
that  the  states  represent  change  between  epochs.  When  pseudolites  go  out  or  come 
back  into  view,  the  number  of  ambiguity  states  changes  from  one  epoch  to  another. 
When  the  double  difference  base  PRN  changes,  the  quantities  that  the  ambiguity 
states  represent  also  change.  To  account  for  these  occurrences,  the  same  functions 
that  were  used  in  the  forward  filter  to  handle  this  in  the  forward  filter  are  adapted 
for  the  optimal  smoothing  algorithm. 

When  the  smoothed  estimate,  x(tj|fj),  is  generated  backward  in  time,  it  is 
formed  from  x(ti+1|£f)  and  H{t~+1),  which  both  correspond  to  the  next  epoch.  If 
pseudolite  4  was  in  view  at  ti+i  but  not  at  U  both  x(ti+i|t/)  and  St(t~+1)  will  have  to 
be  modified  just  like  they  lost  a  pseudolite  clue  to  visibility.  They  are  modified  by 
removing  the  state  estimate  for  pseudolite  4,  along  with  the  corresponding  rows  and 
columns  from  the  covariance  matrix. 

Conversely,  if  pseudolite  4  was  in  view  at  ti  but  not  at  the  states  x(ti+i|t/) 
and  a(t^+1)  will  have  to  add  the  state  estimate  for  this  PRN,  in  addition  to  their 
covariances  P(tj+i|t/)  and  P(t“|„1),  using  the  procedure  shown  in  Section  3.6.2. 
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When  the  base  PRN  is  changed,  the  states  and  covariances  must  also  change 
in  a  similar  manner.  The  states  x(tj+i|t/)  and  x(t~+1)  must  be  altered  to  reflect 
the  correct  double  difference  base  PRN.  The  covariances  are  also  required  to  change 
accordingly. 

3.6.5  Nonlinear  Filtering.  A  number  of  methods  or  techniques  have  been 
proposed  to  deal  with  measurement  nonlinearities  that  were  introduced  in  Section 
2.6.3.  One  was  is  to  simply  increase  the  measurement  variances,  as  done  by  refer¬ 
ences  [9,  32],  A  second  method  is  to  implement  a  full  second  order  filter  based  on  a 
second  order  Taylor  series  approximation  to  the  nonlinear  model.  A  third  method 
has  proposed  to  just  include  the  bias  correction  terms  only  [1,  20].  This  has  been 
shown  to  produce  very  similar  performance  to  the  full  second  order  filter  without 
the  computational  burden  of  the  second  moment  calculations  [20].  The  first  order 
extended  Kalman  filter  update  equations  can  be  modified  to  include  second  order 
terms  for  nonlinear  filtering,  yielding 


where 


A(ti)  =  H(ti)P(tr)Hr(ii)  +  B  m(tr)  +  R{U)  (3.47) 

K  (U)  =  P(tr)HT(ii)A(ti)“1  (3.48) 

x(t+)  =  x(t“)  +  K (ti)  {z (U)  -  h  [x(t“),ti]  -  bm(t~) }  (3.49) 

P (tf)  =  P (tr)  -  K(78)H(it)P(t-)  (3.50) 

(3.51) 


and  the  measurement  bias  correction  term  is  defined  by 


b mkiti  )  =  -tr 


d‘2hk\yi{ti  ),U 
<9x2 


p(f.- 


(3.52) 


with  k  =  1, 2,  •  •  •  ,  m  and  m  represents  the  number  of  measurements. 
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The  differences  between  the  three  modified  nonlinear  filters  are  in  the  deriva¬ 


tion  for  B m{ti  )  where: 


Bias  Correction  Term  Only 

B m{ti  )  =  0 

Modified  Truncated  Second  Order 

B mkl(ti  )  =  —  b mkiti  )b^(tj  ) 

Modified  Gaussian  Second  Order 

B mkiiti  )  =  \tr  <j 

f  d2hk[ic(ti  ),ti]p/ 92/ii[x(C  ),tj]p/,-sl 
9x2  ^  *  )  9x2  v  *  /  J 

1 

The  non-zero  second  partial  derivatives  are  given  as 


d2hk\x.,U 

d2x\ 


x=x(ti  ) 


d'2hk[x,ti 

82X2 


x=x(ti  ) 


d2hk\x,ti 

d2x3 


x=x(ti  ) 


d  2hk\x.,U 

dx\dx2 


d2hk[x,ti 

dx\dxz 


x=x(C  ) 


d2hk[x,ti 

dx2dx3 


x=x(ti  ) 


1  f  ((xJ'-q;i)2+(;/J-a:2)2+(zJ-J3)2)2-(^i-^i)2  1 
^  l  [(xi  —xi^+^yi  —X2)2+(zi  —X3)2]3/2  J 


p  f  ((xfc— a;i)2+(i/*:— X2)2+(zfc— CC3)2)2  — (xfc— a:i)2  1 

^  [  [(xk -xi)2+(yk -X2)2+(zk  —X3)2^2  J 
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^  l  [(xi -xi)2+(yi -x2)2+(zi -X3)2]3/2  J 
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[(a;fe— a;i)2+(j/fc— X2)2+(^fe— a^3)2] 

_ ~(yj -x2)(zj -X3) _ 
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(3.53) 


(3.54) 


(3.55) 


(3.56) 


(3.57) 


(3.58) 
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3.6.6  Weighted  Measurement  Covariance  Matrix.  This  section  describes  a 
novel  method  to  weight  the  measurement  covariance  matrix  R  selectively,  based  on 
the  predicted  tropospheric  delay  for  each  measurement  generated  by  the  tropospheric 
model.  This  allows  the  filter  to  weight  measurements,  relative  to  their  corresponding 
predicted  tropospheric  delay. 

Tropospheric  delay  is  a  function  of  the  atmospheric  parameters,  slant  range, 
and  elevation  angle.  The  atmospheric  parameters  are  constant  for  all  measurements 
in  a  given  simulation  run.  The  longer  the  slant  range  and  the  lower  the  elevation 
angle,  the  larger  the  tropospheric  delay.  If  two  measurements  have  the  same  atmo¬ 
spheric  values,  and  one  measurement  has  a  larger  predicted  delay,  it  is  due  to  either 
a  longer  slant  range  and/or  a  lower  elevation  angle.  For  example,  if  the  tropospheric 
filter  computes  a  tropospheric  delay  of  8  nr  for  prn  1,  and  2  m  for  PRN  2,  the  residual 
tropospheric  error  after  a  model  can  be  expected  to  be  4  times  larger  for  PRN  1  than 
for  PRN  2.  It  would  follow  that,  as  the  predicted  delay  increases,  the  corresponding 
value  in  the  R  matrix  should  also  increase.  The  standard  deviations  of  the  mea¬ 
surements  (i.e.,  the  square  roots  of  the  variance  terms  along  the  diagonal  of  the  R 
matrix)  for  the  phase  is  the  Root  Sum  Square  (RSS)  of  the  standard  deviations  for 
the  tropospheric  delay  error  and  the  non-tropospheric  errors.  If  these  standard  de¬ 
viations  are  0.09  m  and  0.07  m  respectively,  the  RSS  is  0.114  m  or  0.5991  LI  cycles. 
The  variance  is  0.59912  =  .359  and  the  covariance  is  half  the  variance  as  described 
in  Section  3.4.2,  which  results  in  the  baseline  R  matrix  for  the  phase  partition  as 
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This  is  the  approach  used  in  the  baseline  filter  for  the  measurement  covariance  ma¬ 
trix  computation.  When  the  weighted  measurement  covariance  method  is  used,  the 
standard  deviation  for  the  tropospheric  component  is  not  assumed  to  be  0.09  m  for 
every  measurement.  Instead,  the  standard  deviations  for  the  tropospheric  compo¬ 
nents  were  determined  by 


&tropo  =  |-03  x  tropodelay\ 


.161 

.080 

.053 

.213 

.004 


(3.60) 


where  the  absolute  value  of  the  scaled  (0.03)  tropospheric  delay  in  meters  to  result 
in  a  vector  of  standard  deviations.  When  these  values  were  used  in  place  of  a  vector 
of  0.09  m  values  the  R  matrix  which  is  calculated  as  .05  *  (atr0po  *  cr^  ),  is  given  as 
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.213 
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.129 

.328 

.103 
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.214 

.272 

.085 
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.328 

.272  1.384 

.217 
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.103 

.085 

.217 

.136 

The  scale  factor  0.03  was  chosen  because  it  typically  was  equally  likely  to  produce 
values  above  0.09  as  it  was  to  produce  values  below  0.09  for  this  trajectory. 


3.6.7  Tropospheric  Model  Error  States.  The  errors  in  the  tropospheric 
model  include  measurement  errors  in  the  sensors,  atmospheric  errors  due  to  ground 
effects,  and  the  use  of  estimated  positions  of  the  transmitters  and  receivers. 
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The  first  error  source  is  the  set  of  errors  due  to  imprecise  instruments  for 
measuring  atmospheric  pressure,  temperature,  and  relative  humidity.  These  errors 
will  affect  all  measurements  by  roughly  the  same  percentage. 

The  second  error  is  due  to  ground  effects  from  foliage  and  buildings.  The  height 
of  a  typical  test  mission  could  be  2500  meters  above  the  earth’s  surface,  so  ground 
effects  that  only  affect  the  first  25  meters  only  represent  1  percent  of  the  total  signal 
range. 

The  third  error  is  due  to  using  estimated  positions  of  the  transmitters  and 
receivers.  Because  these  errors  are  in  the  centimeter  range,  their  effect  is  almost 
insignificant  to  the  total  error  of  the  tropospheric  model. 

Figure  3.6.7  shows  the  true  tropospheric  delay  in  the  top  subplot  and  the 
residual  tropospheric  error  after  the  model  was  applied  in  the  bottom  subplot.  This 
was  for  a  typical  simulation  run. 
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Figure  3.7  True  Tropospheric  Delay  and  Residual  Tropospheric  Error 
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Although  the  plots  in  Figure  3.7  appear  to  be  the  same  plot,  they  are  not.  Note 
the  scaling  of  the  top  subplot  is  in  meters  and  the  bottom  subplot  is  in  centimeters. 
The  errors  are  highly  correlated  to  the  true  tropospheric  delay.  When  the  errors 
are  expressed  as  a  percentage  of  the  true  delay,  they  are  nearly  constant  and  are 
shown  in  Figure  3.8  for  both  the  mobile  receiver  in  the  top  subplot  and  the  reference 
receiver  in  the  bottom  subplot. 
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Figure  3.8  Residual  Tropospheric  Error  Expressed  as  a  Percent  of  Tropospheric 
Delay 

If  the  residual  tropospheric  errors  after  a  model  is  applied  are  nearly  constant, 
they  can  be  modeled  and  removed.  The  error  percentages  for  both  the  mobile  and 
reference  receiver  were  modeled  as  First  Order  Gauss  Markov  (FOGM)  process.  Two 
states  were  added  to  the  baseline  filter,  one  each  for  the  mobile  and  reference  receiver 
tropospheric  error  percentages.  These  two  states  were  added  after  the  acceleration 
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states  and  before  the  ambiguity  states  and  are  defined  as 


xio  =  (— 1/Tf)xi0  +  wtl  (t) 

xn  =  (-1/Tt)xn  +  wt2(t)  (3.62) 

with  associated  dynamic  driving  noise  processes  given  by 

2  a2 

E  {wtl(t)wtl(t  +  r)}  =  E  {wt2(t)wt2(t  +  t)}  =  — -<5(t)  =  qt5(r)  (3.63) 

J-t 

The  correlation  time,  Tt,  and  variance,  of,  are  determined  based  on  the  anticipated 
error  percentages  and  time  correlations.  The  value  Tt  was  set  to  75  hours  (270000 
seconds)  to  account  for  typical  changes  in  atmospheric  effects,  and  crt  was  set  to  0.03 
m  to  handle  the  ’’worst  case”  error  percentage.  With  these  values  qt  is  calculated  to 
be  1.25  x  10-7m2/sec. 

3.7  Carrier-phase  Ambiguity  Resolution 

The  structure  of  the  carrier-phase  ambiguity  resolution  techniques  are  shown 
in  Figure  3.9.  First  the  Z-Transformation  is  applied  to  the  floating-point  ambiguity 
estimates  and  covariances.  Next,  FASF  generates  the  candidate  ambiguity  sets.  If 
more  than  one  set  is  generated,  a  ratio  test  determines  the  best  ambiguity  set,  based 
on  the  sum  of  square  residuals.  The  inverse  Z-Transformation  is  applied  to  bring  the 
selected  set  of  ambiguities  back  from  the  LAMBDA  domain.  LAMBDA,  FASF,  and 
the  ratio  test  were  described  in  Chapter  2. 


3-38 


Figure  3.9  Ambiguity  Resolution  Algorithm  Description 

3.8  Chapter  Summary 

This  chapter  described  the  truth  model  and  the  error  generation  for  both 
the  measurements  and  measurement  model.  The  floating-point  filter  was  developed 
along  with  some  modifications.  Lastly,  the  structure  of  the  carrier-phase  ambiguity 
resolution  process  was  described. 
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IV.  Results 


4-1  Overview 

This  chapter  presents  simulation  results  and  analysis  of  the  algorithm  devel¬ 
oped  in  this  research.  The  first  section  provides  the  simulation  parameters,  scenario 
descriptions,  and  test  case  definitions  that  are  used  throughout  this  chapter.  The 
next  section  analyzes  a  single  filter  run  to  demonstrate  the  performance.  In  addition 
to  single-run  analysis,  Monte  Carlo  simulations  were  conducted  in  which  an  error 
sensitivity  analysis  and  a  comparison  of  various  filter  enhancements  was  performed. 
Sections  are  also  included  for  optimal  smoothing,  and  alternate  aircraft  trajectories. 
Recall  that  the  primary  purpose  of  this  research  is  to  define  design  tradeoffs  via  sim¬ 
ulations,  and  not  necessarily  to  give  an  absolute  measure  of  the  filter’s  performance. 

4-2  Simulation  Parameters,  Scenario  Descriptions  and  Test  Case  Definitions 

This  section  describes  the  simulated  trajectory,  the  atmospheric  parameters 
used  in  the  truth  model,  along  with  the  scenario  descriptions  and  case  definitions. 
The  statistics  and  criteria  used  to  evaluate  the  filter  performance  are  also  described. 

As  stated  in  Section  3.2.2,  the  trajectory  was  from  actual  C-12  flight  data.  The 
832  second  (14  minute)  section  that  was  used  in  this  research  will  be  referred  to  as 
the  main  flight  trajectory.  This  trajectory  resulted  in  mobile-receiver-to-transmitter 
ranges  of  3  to  32  kilometers.  Other  flight  trajectories  were  also  investigated  and  they 
will  be  described  in  Section  4.6.  For  the  main  flight  trajectory,  the  maximum  range 
from  the  mobile  receiver  to  any  pseudolite  was  48  kilometers,  while  the  maximum 
range  from  the  reference  receiver  to  any  pseudolite  was  32  kilometers.  To  simulate 
pseudolites  coming  into  and  going  out  of  view,  the  maximum  range  allowed  between 
any  pseudolite  and  the  mobile  receiver  was  set  at  32  kilometers  (just  longer  than 
the  maximum  range  to  the  reference  receiver).  Any  measurement  over  this  limit  was 
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considered  out  of  pseudolite  transmission  range.  This  constraint  caused  the  number 
of  visible  pseudolites  at  the  mobile  receiver  to  vary  between  7  and  10  throughout  the 
simulation,  which  was  shown  in  Figure  3.6  from  Section  3.6.1. 

The  atmospheric  parameters  used  for  the  tropospheric  truth  model  are  sum¬ 
marized  Table  4.1.  Atmospheric  pressure  is  the  force  per  unit  area  exerted  against 
a  surface;  by  the  weight  of  the  air  above  that  surface  which  is  sometimes  known  as 
barometric  pressure  [13].  The  average  atmospheric  pressure  at  sea  level  is  1013.25 
millibars  or  29.92  inches  of  mercury.  This  was  the  value  used  for  the  truth  model 
for  every  simulation.  Recall  that  the  filter  model  uses  error-corrupted  atmospheric 
values  to  calculate  tropospheric  delay.  The  value  of  the  temperature  was  chosen  as 
52  degrees  Fahrenheit  because  that  was  approximately  the  average  yearly  tempera¬ 
ture  for  Holloman  AFB,  NM.  Relative  humidity  was  chosen  as  35  percent,  because 
it  is  a  reasonable  value  for  the  desert  climate  at  Holloman. 

Table  4.1  Atmospheric  Parameters  Used  in  Truth  Model 


Atmospheric  pressure 

1013.25  rnbar  (29.92  in) 

Temperature 

284.26  K  (52  F) 

Relative  humidity 

35  percent 

Every  simulation  is  classified  by  a  scenario  description  and  test  case  definition. 
Each  scenario  represents  a  different  objective  which  is  shown  in  Table  4.2.  The  first 
two  scenarios  in  the  table  evaluated  the  filter’s  performance  with  the  baseline  filter 
for  both  the  single  and  widelane  observable.  The  next  three  scenarios  examined  the 
filter’s  sensitivity  to  each  error  source,  with  the  remaining  scenarios  corresponding 
to  each  of  the  filter  enhancements  described  in  Chapter  3. 
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Table  4.2  List  of  Scenario  Descriptions 


Scenario  Identifier 

Scenario  Description 

Base 

Baseline  Filter  with  Single  (LI)  Frequency 

Base  WL 

Baseline  Filter  with  Widelane  Frequency 

PLE 

Impact  of  Pseudolite  Location  Error 

Meas.  Noise 

Impact  of  Measurement  Noise 

Multipath 

Impact  of  Multipath  Error 

Tropo.  Delay 

Impact  of  Tropospheric  Delay 

PR  ave. 

LI  and  L2  Code  Averaging 

Bias  Corr. 

Bias  Correction  Terms 

Smoother 

Optimal  Smoothing 

Weighted  R 

Weighted  Measurement  Covariance  (R)  Matrix 

Tropo.  State 

Tropospheric  Model  Error  States 

All  Enh. 

All  Enhancements  except  Tropo.  State 

The  measurement  errors  that  were  described  in  Chapter  3  represent  the  nomi¬ 
nal  error  case.  In  order  to  test  realistic  magnitudes  of  measurement  errors  fully,  the 
best  and  worst  case  scenarios  were  implemented.  They  represent  one  half  and  twice 
the  nominal  case  for  each  of  the  measurement  error  sources. 

The  Monte  Carlo  simulation  included  evaluation  criteria  for  both  the  floating¬ 
point  results  and  the  ability  to  resolve  the  ambiguities.  The  floating-point  criteria 
included  the  Root  Mean  Square  (RMS)  of  the  three-dimensional  position  error,  the 
RMS  of  the  ambiguity  error,  and  the  RMS  of  the  standard  deviations  for  the  ambi¬ 
guity  state  estimates  (square  root  of  the  state  covariance  matrix  entries). 

The  RMS  position  error  calculation  was  a  three  dimensional  Distance  RMS 
(DRMS)  which  can  also  be  referred  to  as  a  Mean  Radial  Spherical  Error  (MRSE) 
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defined  as 


3D  RMS  = 
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E  (*?  +  yf  +  *?) 
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(4.1) 


where  Xi,  yt,  and  zt  represent  the  position  error,  in  each  axis,  and  m  is  the  number 
of  measurement  epochs.  These  are  temporally  arranged,  versus  ensemble-averaged 
statistics.  The  RMS  statistic  for  the  ambiguity  error  was  defined  as 
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where  Xj(ti )  is  the  ambiguity  error  which  was  summed  for  each  each  ambiguity  per 
epoch  (n)  and  for  each  epoch  (m).  The  RMS  value  for  the  ambiguity  standard 
deviations  were  calculated  in  a  similar  where 
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with  (73{ti)  representing  the  standard  deviation  (square  root  of  the  ambiguity  state 
covariance)  of  a  given  ambiguity.  For  the  Monte  Carlo  simulation,  these  values  were 
evaluated  for  epochs  400  through  832  (final  7  minutes)  to  allow  the  filter  to  converge 
before  calculating  statistics  on  the  performance. 

The  RMS  statistic  was  chosen  to  represent  floating-point  filter  performance  be¬ 
cause  it  required  only  one  value,  versus  two  if  mean  and  standard  deviation  statistics 
were  used.  The  primary  objective  was  the  fixed-integer  carrier-phase  performance. 

In  addition  to  the  floating-point  filter  criteria,  the  percentages  of  correct  fixes 
for  both  simple  rounding  and  the  ambiguity  resolution  techniques  described  in  Chap¬ 
ter  3  are  used.  Because  ambiguity  resolution  cannot  always  resolve  the  ambiguities, 
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the  percentage  of  incorrect  and  unresolved  epochs  are  also  given,  as  shown  in  Table 
4.3. 


Round  %corr 

Percentage  of  correct  fixes  using  simple  rounding 

Amb.  Res.  %corr 

Percentage  of  correct  epochs  with  LAMBDA  and  FASF 

Amb.  Res.  %incorr 

Percentage  of  incorrect  epochs  with  LAMBDA  and  FASF 

Amb.  Res.  %unres 

Percentage  of  unresolved  epochs  with  LAMBDA  and  FASF 

Table  4.3  Ambiguity  Resolution  Evaluation  Criteria 


The  floating-point  ambiguity  state  estimates  and  their  associated  covariances 
were  saved  for  10  equally  spaced  epochs  in  the  second  half  of  the  simulated  test 
run.  The  first  evaluation  criterion  was  to  apply  simple  rounding  of  the  ambiguities, 
and  then  to  calculate  the  percent  of  the  epochs  that  were  correct  (i.e. ,  where  all 
of  the  floating-point  ambiguities  were  within  one-half  of  a  cycle  of  the  true  integer 
ambiguity).  The  rounded  solution  was  considered  correct  only  if  all  ambiguities 
for  that  epoch  were  correct.  The  LAMBDA/FASF  ambiguity  resolution  techniques 
described  in  Section  2.5.8  were  also  applied  to  the  same  floating-point  state  estimates 
and  the  percentage  of  correct  fixes,  percentage  of  incorrect  fixes,  and  percentage  of 
epochs  where  the  ambiguities  were  unresolved,  were  all  computed. 

4-3  Single  Run  Performance 

This  section  evaluates  a  typical  run  of  the  baseline  filter  for  the  single  (LI) 
observable.  The  position  error,  velocity,  acceleration,  measurement  residuals,  and 
ambiguity  errors  are  investigated.  While  this  section  examines  the  filter  for  a  single 
run,  the  following  section  provides  Monte  Carlo  simulation  analysis  in  order  to  pro¬ 
vide  a  complete  picture  of  filter  performance.  Recall  from  Section  3.3  that  the  truth 
model  only  determined  the  true  position  of  the  receivers  along  with  measurement 
errors.  Therefore,  error  plots  are  not  presented  for  the  velocity  and  accelerations. 
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Although  the  primary  goal  of  the  filter  is  to  achieve  a  high  level  of  accuracy  for 
the  ambiguity  estimates,  the  error  in  estimated  receiver  position  can  indicate  how 
well  the  filter  is  performing.  The  X,  Y,  and  Z  position  errors  are  shown  in  Figure 
4.1,  along  with  the  filter-computed  standard  deviations  (square  root  of  the  position 
variances).  Note  that  the  position  errors  are  typically  less  than  the  one  standard 
deviation  (particularly  after  the  initial  transients)  which  may  indicate  that  the  filter 
is  tuned  somewhat  conservatively. 


399849  399988  400126  400265  400403  400542  400680 

08:04  08:06  08:09  08:11  08:13  08:16  08:18 

GPS  Week  Seconds-Local  Time 


Figure  4.1  Position  Errors  and  Filter-Computed  Standard  Deviations  (Dashed 
lines) 

Notice  that  the  time  label  given  in  minutes  are  rounded  to  the  nearest  minute. 
This  label  was  added  to  give  the  reader  a  sense  of  the  relative  time  frame  of  the  test 
run. 

A  3-dimensional  position  error  plot  more  clearly  shows  the  centimeter  level 
accuracy  attained  by  the  filter  in  Figure  4.2.  Note  that  this  still  represents  floating¬ 
point  ambiguities  because  ambiguity  resolution  has  not  taken  place.  Recall  that  a 
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widelane  frequency  aids  the  ambiguity  resolution  process  at  the  expense  of  degraded 
performance  in  the  position  solution.  If  a  widelane  frequency  were  used,  the  perfor¬ 
mance  in  the  position  domain  would  not  achieve  this  same  level.  More  discussion  on 
the  difference  in  accuracies  between  the  widelane  and  single  frequency  is  included  in 
Section  4.5.1. 


GPS  Week  Seconds-Local  Time 


Figure  4.2  3-Dimensional  Position  Error 

The  filter  also  computes  the  velocities  and  accelerations  in  the  ECEF  frame. 
When  the  truth  model  was  generated,  only  the  position  of  the  mobile  receiver  was 
simulated.  It  is  still  beneficial  to  plot  the  velocities  and  accelerations  (Figures  4.3 
and  4.4)  in  order  to  show  when  the  mobile  receiver  experienced  large  accelerations 
(large  for  a  C-12).  Notice  that  the  velocity  and  especially  the  acceleration  plots 
indicate  large  accelerations,  both  in  the  middle  and  at  the  end  of  the  test  run. 
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Residuals  are  an  indicator  of  filter  performance  and  should  be  both  white 
(uncorrelated  in  time)  and  have  a  zero-mean  distribution.  Recall  from  Section  2.2 
that  the  residuals  are  defined  as  the  measurement  minus  the  measurement  prediction. 
The  residuals  for  the  VAp1-2  and  VA01-2  measurements  (for  comparison  purposes) 
are  shown  in  Figure  4.5.  Note  the  scale  of  the  error  for  both  is  given  in  meters, 
with  the  code  plot  having  a  larger  error  range.  The  code  residual  appeared  to  be 
white  for  the  majority  of  the  plot,  with  a  slight  time  correlation  at  the  end  of  the 
simulation  run.  The  phase  residual  typically  displayed  a  much  smaller  magnitude 
than  the  code  residual,  but  showed  more  time  correlation.  The  residuals  are  a 
product  of  measurement  errors  and  dynamics.  The  phase  residuals  have  much  smaller 
measurement  errors,  which  reduces  the  size  of  phase  residuals.  The  dynamics  of 
the  receiver  affect  the  code  and  phase  in  exactly  the  same  manner.  The  phase 
residual  plot  more  clearly  shows  the  effect  of  vehicle  dynamics,  because  unlike  the 
code  measurements,  it  is  not  obscured  by  the  measurement  errors.  These  residual 
errors  are  due  to  the  inability  of  the  full-state  filter  to  predict  future  error  dynamics 
precisely  during  state  propagation. 


4-9 


4 


>  o^V,. 

CD  >?£ 


v  •  '••■* 


.  > 


— w\  ■•; 
’•■■•:  V,' 


-2 1 - 

399849 


399988 

08:06 


400126 

08:09 


400265 

08:11 


400403 

08:13 


400542 

08:16 


400680 

08:18 


GPS  Week  Seconds-Local  Time 


Figure  4.5  Code  and  Phase  Ambiguity  Residuals 


The  ambiguity  error  for  each  of  the  9  ambiguities,  along  with  their  filter- 
computed  1-as  (square  root  of  their  variances),  are  shown  in  Figure  4.6.  Pseudolites 
can  display  much  higher  levels  of  relative  motion  with  their  receivers  than  GPS 
satellite  transmitters  can.  This  will  result  in  the  filter  converging  to  an  ambiguity 
estimate  more  quickly.  This  filter  showed  convergence  with  the  first  5-7  minutes. 
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Figure  4.6  Ambiguity  Error  and  Ambiguity  1-cr  Plots 

4-4  Comparison  of  Code  Versus  Relative  Motion  for  Ambiguity  Convergence 

The  two  factors  that  affect  the  ability  of  the  Elter  to  converge  to  the  correct 
ambiguity  values  are  the  influence  of  the  code  measurements  and  the  relative  motion 
between  transmitters  and  receivers.  This  section  analyzes  how  each  of  these  factors 
affect  filter  convergence. 

The  influence  of  the  relative  motion  was  investigated  by  artificially  increasing 
the  covariance  of  the  code  measurement  to  the  point  that  the  code  measurements 
carry  very  little  weight  in  the  filter.  Recall  from  Section  3.5.2  that  the  covariance 
of  the  double  difference  code  measurement  was  (5m)2.  For  this  test,  that  value  was 
set  to  (1,  000,  000m)2.  This  essentially  removes  the  effect  of  code  measurements  and 
forces  the  filter  to  rely  on  relative  motion  only.  The  results  of  this  test  show  that, 
for  this  simple  case,  the  filter  converged  more  quickly  without  code  than  it  did  with 
code  measurements,  as  shown  in  Figure  4.7.  From  this  plot  is  appears  that  relative 
motion  is  the  primary  factor  in  ambiguity  convergence. 
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Ideally,  bringing  in  code  measurements,  or  additional  measurements,  should 
not  degrade  filter  performance.  Recall  that  the  measurement  error  sources,  especially 
the  multipath,  are  not  zero-mean  or  uncorrelated  in  time.  The  degradation  is  not  a 
simple  mistiming  of  the  code  measurement  covariance,  but  rather,  a  mis-modelling 
of  the  error  source  because  the  filter  is  expecting  a  white,  zero-mean  measurement 
error. 


GPS  Week  Seconds-Local  Time 

Figure  4.7  RMS  Ambiguity  Error  and  Ambiguity  1-er  Plots  for  RCO(fe  =  1012m2 

Notice  that  the  covariances  are  identical  in  the  bottom  subplot.  This  shows 
that  the  code  measurement  variances  are  not  a  strong  influence  on  ambiguity  covari¬ 
ances. 

The  influence  of  code  measurements  was  investigated  by  changing  to  a  station¬ 
ary  trajectory  (i.e. ,  no  motion  in  the  mobile  receiver).  This  change  caused  substantial 
errors  in  the  ambiguity  estimates,  increasing  the  ambiguity  estimation  errors  by  ap¬ 
proximately  a  factor  of  10.  The  errors  in  the  filter-corrupted  ambiguity  estimates 
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for  each  of  the  ten  pseudolites  are  shown  in  Figure  4.8.  This  further  supports  the 
claim  that  relative  motion  is  the  primary  factor  in  ambiguity  convergence. 
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Figure  4.8  Filter  Ambiguity  Estimate  Error  and  Ambiguity  1-cr  Plots  for  Trajec¬ 
tory  without  Motion 

The  analysis  suggests  a  reliance  on  the  relative  motion  between  the  transmit¬ 
ters  and  receivers  for  accurate  ambiguity  estimation.  These  results  suggest  that  a 
pseudolite-based  flight  reference  system  does  not  need  to  depend  much  upon  code 
measurements. 

4-5  Monte  Carlo  Performance 

This  section  first  presents  the  Monte  Carlo  simulation  analysis  of  the  baseline 
filter  for  both  a  single  and  widelane  frequency,  then  an  analysis  of  the  sensitivity  of 
the  filter  to  each  error  source,  and  finally  a  comparison  of  each  filter  enhancement 
described  in  Chapter  3.  The  Monte  Carlo  simulation  conducted  in  this  research 
involved  34  separate  tests  of  100  runs  each. 
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The  baseline  filter  and  filter  enhancement  scenarios  each  had  three  cases:  Best, 
Nominal,  and  Worst.  For  these  scenarios,  all  four  error  sources  were  set  together  as 
the  best,  nominal,  or  worst  case,  in  contrast  to  the  sensitivity  error  analysis  where 
only  one  source  at  a  time  was  adjusted  to  the  best  or  worst  case,  while  the  remaining 
error  sources  were  set  to  the  nominal  values.  This  allowed  the  relative  impact  of  each 
error  source  to  be  evaluated. 

4-5.1  Baseline  Filter.  The  first  two  scenarios  involve  the  baseline  filter  for 
both  the  widelane  and  single  frequency  filters.  Each  filter  was  evaluated  against  the 
best,  nominal,  and  worst  error  cases,  which  are  presented  in  Table  4.3. 

Table  4.4  Widelane  versus  Single  Frequency  in  Baseline  Filter 


Test 

# 

Scenario 

Error 

Case 

3-D 

RMS 

VA  Nerr 

RMS 

°VA  N 

RMS 

Round 

%corr 

Amb. 

Res. 

%corr 

Amb. 

Res. 

%incorr 

Amb. 

Res. 

%unres 

1 

WL 

Best 

.063 

.031 

.046 

100 

100 

0 

0 

2 

Single 

Best 

.033 

.051 

.083 

99.4 

99.3 

.6 

.1 

3 

WL 

Nom. 

.123 

.065 

.046 

100 

100 

0 

0 

4 

Single 

Nom. 

.054 

.101 

.083 

99.4 

96.3 

.5 

3.2 

5 

WL 

Worst 

.241 

.126 

.046 

99.2 

94.8 

.3 

4.9 

6 

Single 

Worst 

.099 

.225 

.083 

75.2 

59.2 

.4 

40.4 

The  use  of  a  widelane  frequency  improved  ambiguity  resolution  ability,  but  at 
the  expense  of  a  larger  position  error.  This  is  reasonable  because,  as  shown  in  Section 
2.5.7,  widelaning  increases  the  effect  of  error  sources  that  are  not  frequency  correlated 
(multipath  and  measurement  noise),  but  does  not  change  the  effect  of  error  sources 
that  are  frequency  correlated  (pseudolite  position  errors  and  tropospheric  delay). 
The  position  errors  and  tropospheric  delay  are  reduced  when  expressed  in  cycles, 
which  reduces  the  ambiguity  search  space.  This  is  clearly  shown  by  comparing  the 
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single  frequency  to  widelane  results  for  each  error  case.  The  widelane  tests  resulted 
in  larger  position  domain  errors,  but  smaller  ambiguity  errors.  The  remaining  Monte 
Carlo  simulation  use  single  frequency  (LI)  phase  measurements  in  order  to  show  the 
contribution  of  each  error  source  better,  and  the  performance  improvement  of  each 
filter  enhancement. 

The  floating-point  filter  was  tuned  for  the  nominal  error  case  for  both  the  single 
and  widelane  observables.  The  agreement  between  the  RMS  ambiguity  error  (.101 
from  test  4)  and  RMS  standard  deviation  (.083  from  test  4)  indicate  how  well  the 
filter  is  tuned.  The  results  also  show  that  the  filter  is  too  conservative  for  the  best 
case,  and  is  overestimating  its  performance  for  the  worst  case  error. 

Tests  1,  2,  5,  and  6  represent  a  mistimed  filter  and  therefore  can  be  used  to 
determine  how  robust  the  filter  is  to  mistiming.  These  tests  indicate  that  the  filter 
is  fairly  robust  to  mistiming.  The  percentage  of  correct  fixes  do  decrease  going  from 
the  nominal  error  case  to  the  worst  case,  but  the  percent  of  incorrect  fixes  does  not 
significantly  change.  In  fact,  it  actually  decreased  for  the  single  frequency  case.  The 
biggest  change  is  the  percentage  of  unresolved  cases,  which  could  be  used  to  indicate 
a  mistiming. 

To  investigate  the  capabilities  of  this  algorithm  fully,  the  Liter  was  re-tuned 
for  both  the  best  and  worst  case  single  frequency  filter.  The  re-tuning  of  the  Liter 
involved  increasing  the  measurement  covariance  values  by  a  factor  of  2.5  for  the 
worst  case  and  decreasing  the  values  by  a  factor  of  2.5  for  the  best  case.  The  results 
are  shown  in  Table  4.5,  along  with  Tests  2  and  6  for  comparison.  The  ambiguity 
resolution  process  uses  a  value  of  ”k”  to  determine  how  many  standard  deviations 
dehne  the  search  area.  The  value  of  k  for  the  worst  case  was  set  to  5  instead  of 
10  for  this  test  only.  This  was  due  to  a  large  number  of  unresolved  Lxes.  The 
best  case  re-tuned  Liter  (Test  7)  had  slightly  larger  errors  in  position  and  Loating- 
point  ambiguities,  but  an  improved  ability  to  resolve  ambiguities.  The  worst  case 
re-tuned  Liter  (Test  8)  also  showed  an  slightly  better  performance  in  the  ability  to 
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resolve  ambiguities  in  addition  to  slightly  improving  the  floating-point  position  and 
ambiguity  estimates.  To  summarize,  the  filter  will  perform  best  when  it  is  properly 
tuned,  but  this  filter  can  tolerate  mistiming  to  provide  adequate  results. 


Test 

# 

Scenario 

Error 

Case 

3-D 

RMS 

VA  Nerr 

RMS 

OVA7V 

RMS 

Round 

%corr 

Amb. 

Res. 

%corr 

Amb. 

Res. 

%incorr 

Amb. 

Res. 

%unres 

2 

Single 

Best 

.033 

.051 

.083 

99.4 

99.3 

.6 

.1 

7 

Retuned 

Best 

.036 

.062 

.054 

100 

100 

0 

0 

6 

Single 

Worst 

.099 

.225 

.083 

75.2 

59.2 

.4 

40.4 

8 

Retuned 

Worst 

.097 

.222 

.167 

79.2 

63.8 

1 

35.2 

Table  4.5  Baseline  Single  Frequency  Retuned 


The  baseline  filter  test  showed  a  widelane  frequency  outperforms  the  single 
frequency  in  the  ability  to  resolve  ambiguities,  but  a  single  frequency  is  better  in 
the  position  domain.  The  filter  will  perform  best  when  it  is  properly  tuned,  but  this 
filter  can  tolerate  mistiming  to  provide  adequate  results.  If  the  filter  is  producing  a 
large  number  of  unresolved  epochs,  that  might  indicate  the  measurement  errors  are 
larger  than  anticipated,  which  results  in  the  filter  overestimating  its  performance.  If 
the  filter  is  mistimed,  it  is  better  to  underestimate  than  to  overestimate  its  ability. 

4-5.2  Error  Sensitivity  Analysis.  This  section  examines  the  error  sensitiv¬ 
ity  of  each  error  source.  The  sensitivity  analysis  was  conducted  by  comparing  the 
nominal  case  to  the  best  and  worst  case  of  selected  error  sources.  The  error  source 
of  interest  is  set  to  either  the  best  or  worst  case  error,  while  the  remaining  errors  are 
set  at  the  nominal  values.  The  results  are  compared  to  the  nominal  error  case  (Test 
4).  The  results  for  the  best  case  are  shown  in  Table  4.6. 
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Test 

# 

Scenario 

Error 

Case 

3-D 

RMS 

VA  Nerr 

RMS 

ffVAF 

RMS 

Round 

%corr 

Amb. 

Res. 

%corr 

Amb. 

Res. 

%incorr 

Amb. 

Res. 

%unres 

4 

Single 

Nom. 

.054 

.101 

.083 

99.4 

96.3 

.5 

3.2 

9 

PLE 

Best 

.041 

.075 

.083 

99.4 

99 

.6 

.4 

10 

Meas.  Noise 

Best 

.051 

.100 

.083 

99.4 

97.4 

.5 

2.1 

11 

Multipath 

Best 

.053 

.103 

.083 

99.5 

97.1 

0 

2.9 

12 

Tropo.  Delay 

Best 

.051 

.105 

.083 

100 

98 

0 

2 

Table  4.6  Best  Case  Error  Sensitivity 


The  errors  sources  that  clearly  stand  out  as  the  most  sensitive  are  the  pseu- 
dolite  position  errors  (Test  9)  along  with  the  tropospheric  delay  (Test  12).  This 
is  reasonable  because  these  errors  are  larger  in  magnitude  than  the  multipath  and 
measurement  noise,  in  addition  to  being  more  time  correlated.  When  the  pseudolite 
positions  (Test  9)  were  reduced,  the  percentage  of  correct  fixes  was  99,  which  was 
2.7  percent  higher  than  the  base  case.  The  pseudolite  position  errors  also  showed 
the  only  significant  improvement  in  the  RMS  position  and  RMS  ambiguity  errors. 
The  best  case  tropospheric  delay  (Test  12)  also  showed  an  improvement  in  ambi¬ 
guity  resolution.  The  best  case  for  tropospheric  delay  also  resulted  in  100  percent 
correct  for  the  ambiguity  rounding,  which  means  that  every  floating-point  ambigu¬ 
ity  estimate  was  within  half  a  cycle  of  the  correct  ambiguity.  Both  the  multipath 
and  tropospheric  delay  scenario  achieved  zero  percent  incorrect  while  increasing  the 
percent  correct  over  the  baseline  filter. 

A  sensitivity  analysis  also  examined  the  effect  of  changing  one  error  source  at 
a  time  to  the  worst  case  expected  error.  The  worst  case  error  sensitivity  results  are 
shown  in  Table  4.6.  Each  test  (13-16)  should  show  a  degraded  performance  from  the 
base  (Test  4),  with  the  magnitude  of  degradation  indicating  the  relative  sensitivity. 
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Test 

# 

Scenario 

Error 

Case 

3-D 

RMS 

VA  Nerr 

RMS 

O'VAN 

RMS 

Round 

%corr 

Amb. 

Res. 

%corr 

Amb. 

Res. 

%incorr 

Amb. 

Res. 

%unres 

4 

Single 

Nom. 

.054 

.101 

.083 

99.4 

96.3 

.5 

3.2 

13 

PLE 

Worst 

.083 

.186 

.083 

89.6 

74.9 

.4 

24.7 

14 

Meas.  Noise 

Worst 

.061 

.119 

.083 

99.4 

94 

.4 

5.6 

15 

Multipath 

Worst 

.055 

.115 

.083 

99.4 

95.3 

.6 

4.1 

16 

Tropo.  Delay 

Worst 

.066 

.138 

.083 

95.3 

91.9 

0 

8.1 

Table  4.7  Worst  Case  Error  Sensitivity 


The  pseudolite  position  error  (Test  13)  and  the  tropospheric  delay  (Test  16) 
again  showed  the  most  significant  sensitivity.  Only  the  pseudolite  position  error 
displayed  a  much  larger  change  in  the  RMS  position  and  floating-point  ambiguity 
errors  in  comparison  to  the  other  scenarios. 

The  pseudolite  position  errors  and  the  tropospheric  delay  were  shown  from 
Chapter  3  to  have  larger  magnitudes  and  stronger  time  correlations  than  the  mea¬ 
surement  noise  and  multipath.  As  a  result,  the  filter  is  more  sensitive  to  these  errors 
when  examining  the  best  and  worst  cases  separately  for  each  error  source.  Only  the 
position  errors  impacted  RMS  position  and  RMS  ambiguity  errors  significantly. 

4-5.3  Filter  Enhancements.  The  five  filter  enhancements  that  were  de¬ 
veloped  in  this  research  included  code  averaging,  bias  correction  terms,  optimal 
smoothing,  measurement  covariance  weighting,  and  tropospheric  model  error  states. 
This  section  evaluates  the  filter  enhancements  against  the  best,  nominal,  and  worst 
case  errors  with  a  single  (LI)  frequency.  Tests  2,  4,  and  6  are  used  as  the  baseline 
filter  for  comparison.  Note  that  the  filter  tuning  is  the  same  for  best,  nominal,  and 
worst  case  errors  (i.e.  the  re-tuned  filter  was  not  used). 
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The  tropospheric  model  error  states  will  be  analyzed  for  a  single  run  before 
describing  the  Monte  Carlo  results.  Recall  from  Chapter  3  that  there  were  two 
implementations  developed — two  separate  states  for  the  mobile  and  reference  receiver 
error  percentages,  and  a  single  state  for  the  combination  error  percentage.  When 
the  only  measurement  error  was  the  tropospheric  delay,  both  filters  performed  well, 
which  can  be  seen  in  Figure  4.9. 


GPS  Week  Seconds-Local  Time 


Figure  4.9  Mobile  and  Reference  Rcvr.  Estimated  Tropospheric  Model  Error  Per¬ 
centages  for  Tropospheric  Delay  Only 

Notice  that  the  state  estimates  are  plotted  as  points,  but  appear  as  a  thicker 
line  than  the  plot  of  the  correct  value.  The  next  two  figures  were  plotted  in  the  same 
manner. 

When  the  other  measurement  errors  are  added,  the  filter  is  still  able  to  estimate 
the  mobile  percentage  successfully,  but  not  the  reference  percentage,  as  shown  in 
Figure  4.10. 
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Figure  4.10  Mobile  and  Reference  Rcvr.  Tropospheric  Model  Error  Percentages 
with  all  Measurement  Errors 

At  first  it  was  thought  that  the  reference  error  percentage  was  un-observable, 
dne  to  lack  of  relative  motion  between  the  reference  receiver  and  psendolites.  How¬ 
ever,  if  this  were  true,  the  filter  would  not  have  been  able  to  estimate  the  percentage 
when  only  the  tropospheric  delay  was  present.  Instead,  the  lack  of  relative  motion 
for  the  reference  receiver  means  that  the  psendolite  position  errors  cause  biases  in  the 
measurements.  This  obscures  the  filter  from  correctly  discerning  the  tropospheric 
model  scale  factor  error,  which  looks  like  pseudolite  location  biases.  This  motivated 
the  single  tropospheric  state  as  previously  discussed  in  Chapter  3.  When  this  im¬ 
plementation  was  used,  the  filter-estimated  percentages  typically  fell  between  the 
correct  mobile  and  reference  receiver  percentages,  as  shown  in  Figure  4.11.  The 
single  tropospheric  state  implementation  was  used  for  the  Monte  Carlo  simulations. 
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Figure  4.11  Combined  Tropospheric  Model  Error  Percentage  with  all  Measure¬ 
ment  Errors 

The  correct  mobile  percentage  was  -1  percent  and  the  correct  reference  per¬ 
centage  was  -1.5  percent.  This  plot  shows  the  filter  converging  to  a  value  close  to 
the  correct  mobile  percentage. 

Table  4.8  shows  the  results  for  the  best  error  case.  Both  of  the  enhancements 
that  targeted  the  tropospheric  delay  improved  the  ambiguity  resolution  process.  The 
filter’s  performance  was  very  high  for  the  best  case  and  there  was  relatively  little 
room  for  improvement.  Note  that  the  optimal  smoother  actually  resulted  in  slightly 
worse  results  than  the  baseline  filter.  After  the  filter  has  converged,  typically  there 
is  little,  if  any,  benefit  from  smoothing  due  to  small  dynamics  driving  noises.  The 
results  in  the  table  indicate  the  same  10  points  in  time  at  which  the  smoother  did 
not  outperform  the  filter.  If  another  10  points  in  time  are  selected,  the  smoother 
could  have  outperformed  the  filter.  Basically,  after  a  filter  has  converged  and  has  a 
really  good  dynamics  model,  a  smoother  cannot  be  expected  to  outperform  a  filter  in 
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a  consistent,  significant  manner.  Only  the  tropospheric  model  error  states  increased 
the  accuracy  of  the  filter  in  the  position  domain. 


Test 

# 

Scenario 

Error 

Case 

3-D 

RMS 

VA  Nerr 

RMS 

°VA7V 

RMS 

Round 

%corr 

Amb. 

Res. 

%corr 

Amb. 

Res. 

%incorr 

Amb. 

Res. 

%unres 

2 

Single 

Best 

.033 

.051 

.083 

99.4 

99.3 

.6 

.1 

17 

PR  Ave. 

Best 

.033 

.049 

.083 

99.4 

98.3 

.7 

1 

18 

Bias  Corr. 

Best 

.037 

.052 

.083 

99.4 

99.3 

.6 

.1 

19 

Smoother 

Best 

.041 

.049 

.083 

99.4 

98.4 

.8 

.8 

20 

Weighted  R 

Best 

.033 

.050 

.069 

100 

100 

0 

0 

21 

Tropo.  State 

Best 

.026 

.051 

.089 

100 

99.8 

0 

.2 

Table  4.8  Best  Case  Filter  Enhancement 


Table  4.9  shows  the  results  for  the  nominal  error  case.  The  code  averaging, 
weighted  measurement  covariance,  and  the  tropospheric  model  error  states  could 
have  resolved  100  percent  of  the  ambiguities  correctly  with  simple  rounding.  The 
smoother  offered  the  most  improvement  for  when  comparing  the  percentage  of  cor¬ 
rect  ambiguities  fixes,  but  also  slightly  increased  the  percentage  of  incorrect  Exes. 
Recall  that  the  smoother  actually  decreased  the  performance  for  the  best  case.  The 
code  averaging,  weighted  measurement  covariance,  and  the  tropospheric  model  error 
states  offered  a  more  modest  increase  of  correct  ambiguities,  but  with  zero  percent 
of  incorrect  fixes.  Again,  the  bias  correction  terms  did  not  improve  the  filter’s  per¬ 
formance. 
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Test 

# 

Scenario 

Error 

Case 

3-D 

RMS 

VA  Nerr 

RMS 

O'VAN 

RMS 

Round 

%corr 

Amb. 

res. 

%corr 

Amb. 

res. 

%incorr 

Amb. 

res. 

%unres 

4 

Single 

Nom. 

.054 

.101 

.083 

99.4 

96.3 

.5 

3.2 

22 

PR  Ave. 

Nom. 

.056 

.111 

.083 

100 

96.6 

0 

3.4 

23 

Bias  Corr. 

Nom. 

.057 

.110 

.083 

99.3 

96.3 

.5 

3.2 

24 

Smoother 

Nom. 

.057 

.107 

.083 

99.4 

99.4 

.6 

0 

25 

Weighted  R 

Nom. 

.053 

.105 

.069 

100 

97.9 

0 

2.1 

26 

Tropo.  State 

Nom. 

.041 

.106 

.089 

100 

97.3 

0 

2.7 

Tabic  4.9  Nominal  Case  Filter  Enhancement 


Table  4.10  shows  the  filter  enhancements  results  for  the  worst  case  errors. 
The  optimal  smoother  (Test  29)  and  the  weighted  measurement  covariance  matrix 
(Test  30)  resulted  in  the  highest  percentage  of  correct  fixes,  but  it  also  increased 
the  percent  of  incorrect  fixes.  The  increase  in  percentage  of  correct  fixes  was  27.9 
for  optimal  smoothing  and  14.7  for  the  weighted  measurement  covariance  matrix. 
Again,  code  averaging  (Test  27)  and  the  bias  correction  terms  (Test  28)  showed 
little  improvement  over  the  baseline  filter.  The  tropospheric  model  error  state  was 
the  only  case  to  show  improvement  in  the  position  accuracy.  It  also  had  the  best 
performance  for  simple  rounding,  although  the  ambiguity  resolution  improvement 
was  slight. 
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Test 

# 

Scenario 

Error 

Case 

3-D 

RMS 

VA  Nerr 

RMS 

OVA7V 

RMS 

Round 

%corr 

Amb. 

Res. 

%corr 

Amb. 

Res. 

%incorr 

Amb. 

Res. 

%unres 

6 

Single 

Worst 

.099 

.225 

.083 

75.2 

59.2 

.4 

40.4 

27 

PR  Ave. 

Worst 

.101 

.227 

.083 

76.7 

59.8 

.2 

40 

28 

Bias  Corr. 

Worst 

.101 

.226 

.083 

75.7 

59.2 

.5 

40.3 

29 

Smoother 

Worst 

.095 

.221 

.083 

78.9 

86.9 

2.9 

10.2 

30 

Weighted  R 

Worst 

.098 

.220 

.069 

79.1 

73.9 

3.3 

22.8 

31 

Tropo.  State 

Worst 

.078 

.217 

.089 

82.2 

61.7 

.3 

38 

Tabic  4.10  Worst  Case  Filter  Enhancement 


When  compiling  the  results  for  all  three  error  cases,  the  optimal  smoothing 
and  the  weighted  measurement  covariance  matrix  made  the  largest  difference  in  the 
ability  to  resolve  ambiguities.  Although  they  increased  the  percent  incorrect  in  the 
worst  error  case,  that  could  be  improved  by  proper  tuning  or  by  adjusting  the  ra¬ 
tio  test  criteria  in  the  ambiguity  set  selection  process.  Even  though  code  averaging 
only  slightly  increased  the  performance,  it  is  still  worth  implementing  in  an  oper¬ 
ational  system  in  which  two  signals  are  available.  It  requires  little  computational 
time  and  provides  a  modest  increase  in  the  accuracy  of  the  code  measurements.  The 
tropospheric  model  error  state  method  developed  in  this  research  did  improve  ambi¬ 
guity  resolution,  but  not  to  the  degree  that  smoothing  or  the  weighted  measurement 
covariance  matrix  was  able  to  accomplish.  As  stated  previously,  the  tropospheric 
model  error  state  method  was  the  only  enhancement  to  reduce  errors  in  the  position 
domain  significantly,  and  it  had  the  highest  percentage  of  correct  ambiguities  using 
simple  rounding.  This  suggests  that  this  enhancement  could  possibly  outperform 
the  other  enhancements  with  better  tuning  parameters  in  the  ambiguity  resolution 
process.  The  long  measurement  ranges  for  this  trajectory  did  not  contain  harsh  mea¬ 
surement  model  nonlinearities.  This  would  explain  why  the  bias  correction  terms 
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did  not  significantly  aid  the  ambiguity  resolution  process.  The  next  section  contains 
alternate  trajectories  that  contain  higher  measurement  model  nonlinearities,  with  a 
comparison  of  nonlinear  filtering  techniques  to  address  them. 

The  previous  tests  show  the  benefit  of  each  filter  enhancement,  but  do  not 
show  the  level  of  performance  when  all  the  enhancements  are  used  together.  Ta¬ 
ble  4.11  shows  the  resulting  level  of  improvement  when  all  the  enhancements  are 
used,  except  for  the  tropospheric  model  error  states.  The  weighted  measurement 
covariance  matrix  method  cannot  be  applied  at  the  same  time  the  filter  is  using  the 
tropospheric  model  error  states,  and  the  weighted  measurement  covariance  matrix 
was  chosen  because  it  aided  the  ambiguity  resolution  process  to  a  larger  degree.  The 
table  used  the  worst  error  case  because  it  provided  the  case  with  the  largest  room 
for  enhancement. 


Test 

# 

Scenario 

Error 

Case 

3-D 

RMS 

VA  Nerr 

RMS 

°VA  N 

RMS 

Round 

%corr 

Amb. 

Res. 

%corr 

Amb. 

Res. 

%incorr 

Amb. 

Res. 

%unres 

6 

Single 

Worst 

.099 

.225 

.083 

75.2 

59.2 

.4 

40.4 

32 

All  Enh. 

Worst 

.093 

.216 

.069 

80.8 

76.3 

2.4 

21.3 

Table  4.11  All  Enhancement  Test  for  Worst  Case 


It  is  surprising  to  note  that  the  baseline  filter  with  only  a  smoother  performed 
better  than  Test  32,  where  other  enhancements  were  also  used.  Further  tests  were 
conducted  to  explain  this  phenomenon  by  adjusting  the  scaling  factor  in  the  weighted 
measurement  covariance  (see  Equation  3.60).  When  the  scaling  factor,  which  was 
previously  .03,  was  increased  to  .05,  the  filter  (with  both  the  smoother  and  weighted 
measurement  covariance  matrix)  correctly  resolved  84.5  percent  of  the  ambiguities 
(with  an  incorrect  percentage  of  2.6  percent).  When  the  scaling  factor  increased,  the 
RMS  ambiguity  covariance  also  increased  from  .069  cycles  to  .079  cycles,  which  was 
much  closer  to  the  baseline  filter  RMS  ambiguities  covariance  of  .083  cycles.  This 
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shows  that  the  weighted  measurement  covariance  matrix  is  sensitive  to  the  scaling 
factor. 

This  test  alone  would  suggest  that  the  scaling  factor  should  be  increased  for  all 
the  tests  using  the  weighted  measurement  covariance.  Recall  that  the  worst  error  case 
effectively  has  filter  mistiming,  and  the  lower  ambiguity  covariance  from  the  weighted 
measurement  covariance  matrix  enhancement  would  represent  a  higher  degree  of 
mistiming  than  the  baseline  filter.  This  mistiming  explains  why  the  smoother  had 
better  performance  without  a  weighted  measurement  covariance  for  the  worst  error 
case.  In  order  to  prove  this,  two  additional  tests  were  conducted  with  exactly  the 
same  parameters  as  Test  25,  except  with  different  scaling  factors  in  the  weighted 
measurement  covariance  formulation.  The  results  arc  listed  in  Table  4.12  for  scaling 
factors  of  .05  and  .06,  along  with  the  original  scaling  factor  of  .03  (Test  25).  It  is  clear 
that  increasing  the  scaling  factor  for  the  nominal  error  case  degrades  performance. 
These  results  support  the  claim  that  the  lack  of  performance  in  the  worst  error  case 
was  due,  in  large  part,  to  filter  mistiming.  It  is  important  to  note  that  the  weighted 
measurement  covariance  greatly  increased  accuracy  of  ambiguity  resolution  in  all 
three  error  cases. 


Test 

# 

Scenario 

Error 

Case 

3-D 

RMS 

VA  Nerr 

RMS 

O'VAJV 

RMS 

Round 

%corr 

Amb. 

Res. 

%corr 

Amb. 

Res. 

%incorr 

Amb. 

Res. 

%unres 

25 

Weighted  .03 

Nom. 

.053 

.105 

.069 

100 

97.9 

0 

2.1 

33 

Weighted  .05 

Nom. 

.057 

.108 

.079 

99.4 

87.7 

.5 

11.8 

34 

Weighted  .06 

Nom. 

.059 

.110 

.084 

99.2 

86.0 

.6 

13.4 

Table  4.12  Weighted  Measurement  Covariance  Matrix  Scaling  Factor  Comparison 


The  value  of  r  for  the  tropospheric  model  error  states  was  set  to  75  hours  for 
the  previous  cases,  which  was  discussed  in  Chapter  3.  This  tuning  parameter  was 
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also  set  to  5  minutes  and  1  sec  to  evaluate  its  impact,  with  results  shown  in  Table 
4.13. 


Test 

# 

Scenario 

Error 

Case 

3-D 

RMS 

V  NN err 

RMS 

O'VAN 

RMS 

Round 

%corr 

Amb. 

Res. 

%corr 

Amb. 

Res. 

%incorr 

Amb. 

Res. 

%unres 

26 

Tropo  75  h 

Nom. 

.041 

.106 

.089 

100 

97.3 

0 

2.7 

33 

Tropo  5  m 

Nom. 

.045 

.108 

.094 

100 

95.7 

0 

4.3 

34 

Tropo  1 s 

Nom. 

.048 

.109 

.088 

99.4 

96.1 

.6 

3.3 

Table  4.13  Tropospheric  Model  Error  States  Time  Constant  Comparison 


When  the  time  constant  decreased,  the  values  in  the  tropospheric  model  error 
state  (xio)  appeared  to  be  more  uncorrelated  in  time.  As  the  time  constant  decreases, 
the  filter  is  effectively  allowing  the  state  estimate,  which  is  absorbing  some  of  the 
multipath  and  pseudolite  location  errors,  to  vary  more  with  time.  This  is  important 
because  the  errors  themselves  vary  with  time. 

4-6  Alternate  Trajectories  with  Nonlinear  Filter  Comparisons 

This  section  compares  three  different  trajectories  and  analyzes  the  impact  of 
measurement  model  nonlinearities  for  each  one.  The  evaluation  for  each  trajectory 
involves  investigation  of  nonlinear  filtering  techniques  to  include  1)  EKF  with  bias 
correction  terms,  2)  a  modified  truncated  second  order  filter,  and  3)  a  modihed 
Gaussian  second  order  filter  [20]. 

4-6.1  Landing  Scenario.  The  previous  simulations  used  a  flight  trajectory 
that  stayed  approximately  3  kilometers  above  the  pseudolite  network.  The  slant 
ranges  for  this  trajectory  were  3  to  32  kilometers.  At  these  ranges,  the  bias  correction 
terms  did  not  show  significant  improvement  over  the  EKF  in  terms  of  ambiguity 
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estimation  and  ambiguity  resolution  performance.  This  section  investigates  a  landing 
scenario  and  the  impact  of  measurement  model  nonlinearities. 

The  trajectory  is  actually  the  takeoff  of  a  C-12,  which  was  reversed  to  emulate 
a  landing  scenario.  Figure  4.12  shows  the  3-D  trajectory  (with  projections  onto  each 
axis),  while  Figure  4.13  shows  the  same  pseudolite  and  reference  receiver  locations 
with  the  ground  projection  of  this  new  trajectory.  The  maximum  range  for  the 
pseudolite  signals  was  again  set  to  32  kilometers.  This  trajectory  resulted  in  9-10 
pseudolites  in  view.  It  is  important  to  note  that,  in  this  trajectory,  the  shortest 
range  to  the  mobile  receiver  was  2.5  kilometers. 


Figure  4.12  Landing  Scenario  Trajectory  Plot 
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Figure  4.13  Landing  Scenario  Pseudolite  and  Reference  Receiver  Locations  with 
Trajectory  Ground  Projection 

The  results  shown  in  Table  4.13  show  that  only  a  slight  advantage  is  gained 
with  bias  correction  terms  or  second  order  filtering,  with  the  majority  of  the  benefit 
coming  from  the  bias  correction  terms.  Note  that  the  bias  correction  terms  alone 
did  better  than  either  second  order  filter. 

Table  4.14  Widelane  versus  Single  Frequency  in  Baseline  Filter 


Test 

# 

Scenario 

Error 

Case 

3-D 

RMS 

V  ANerr 

RMS 

&VAN 

RMS 

Round 

%corr 

Amb. 

Res. 

%corr 

Amb. 

Res. 

%incorr 

Amb. 

Res. 

%unres 

35 

EKF 

Nom. 

.298 

.129 

.084 

98.9 

99.7 

0 

.3 

36 

Bias 

Nom. 

.430 

.129 

.084 

99 

99.8 

0 

.2 

37 

Trun. 

Nom. 

.438 

.130 

.083 

98.9 

99.8 

0 

.2 

38 

Gaus. 

Nom. 

.443 

.130 

.084 

98.7 

99.8 

0 

.2 
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This  landing  trajectory  also  did  not  fully  challenge  the  baseline  EKF,  because 
the  shortest  ranges  were  still  in  excess  of  2.5  kilometers.  These  tests  are  still  of 
general  value,  however,  because  they  demonstrate  that  the  algorithm  performs  well 
with  a  different  trajectory. 

4-6.2  Alternate  Landing  Scenario.  Further  test  were  conducted  to  attempt 
to  challenge  the  filter  with  severe  measurement  model  nonlinearities.  The  simulated 
flight  trajectory  was  shifted  in  all  three  directions  to  bring  the  end  of  the  flight  within 
100  meters  of  PRN  8  for  the  last  226  seconds  of  the  832  second  mission  (with  the 
last  192  seconds  at  exactly  73  meters  due  to  a  stationary  trajectory).  This  alternate 
landing  scenario  is  shown  in  Figure  4.14. 
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Figure  4.14  Alternate  Landing  Scenario  Pseudolite  and  Reference  Receiver  Loca¬ 
tions  with  Trajectory  Ground  Projection 

The  results  of  one  simulation  run  for  this  trajectory  are  shown  in  Figures  4.15 
and  4.16  for  the  EKF  and  EKF  with  bias  correction  terms.  The  EKF-only  case 
showed  a  slight  increase  in  ambiguity  error  for  the  end  of  the  run.  The  addition  of 
the  bias  correction  terms  did  not  significantly  improve  the  performance  of  the  filter. 
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The  modified  Gaussian  second  order  filter  was  very  similar  to  the  EKF  with  bias 
correction  terms. 


Figure  4.15 


Figure  4.16 
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Single  Simulation  Run  for  Alternate  Landing  Trajectory  with  an  EKF 
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Single  Simulation  Run  for  Alternate  Landing  Trajectory  with  an  EKF 
with  Bias  Correction  Terms 
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Note  that  one  of  the  ambiguities  in  Figure  4.16  is  worse  than  it  was  in  Figure 
4.15.  This  suggests  that  the  addition  of  the  bias  correction  terms  actually  degraded 
the  performance  at  the  start  of  the  run. 

4-6.3  Take  Off  Scenario.  The  alternate  landing  trajectory  was  inverted 
to  simulated  an  aircraft  taking  off  and  experiencing  the  harsh  nonlinearities  right 
from  the  start  of  the  simulation.  This  forces  the  filter  to  deal  with  the  nonlinearities 
before  the  filter  has  converged.  The  EKF  in  this  scenario  did  require  more  time  to 
converge  and  in  that  time  it  experienced  the  large  ambiguity  errors  shown  in  Figure 
4.17. 
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Figure  4.17  Single  Simulation  Run  for  Takeoff  Trajectory  with  an  EKF  (Initialized 
73  meters  away  from  a  Pseudolite) 

When  the  bias  correction  terms  were  added,  the  filter  diverged  very  quickly 
(within  15  seconds).  The  bias  correction  factors  are  a  function  of  the  covariance 
matrix.  In  the  baseline  filter,  the  position  covariance  values  used  a  standard  deviation 
of  100  meters.  The  filter  diverged  because  the  position  standard  deviation  was 
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greater  than  the  actual  distance  between  the  transmitter  and  receiver,  which  caused 
problems  with  the  line  of  sight  vectors. 

Two  fixes  for  this  problem  were  attempted.  The  first  was  to  delay  the  use  of  the 
bias  correction  terms  until  after  20  seconds  had  elapsed,  to  allow  the  covariance  values 
to  decrease.  This  was  successful,  and  it  produced  results  nearly  identical  to  the  EKF 
case  without  bias  correction  terms.  Lowering  the  covariance  to  a  standard  deviation 
of  4  meters  was  the  second  attempted  solution.  This  also  solved  the  divergence 
problem. 

For  the  application  to  require  any  second  order  nonlinear  filtering  (full-state 
or  just  bias  correction  terms),  the  position  error  must  be  significant  in  relation  to 
the  range  between  the  transmitter  and  receiver.  Second-order  filtering  must  be 
applied  with  care,  because  it  is  sensitive  to  the  covariance  initialization.  When  the 
bias  correction  terms  are  used  with  large  covariance  values,  divergence  can  occur. 
The  modified  truncated  and  modified  Gaussian  second-order  filters  did  not  provide 
sufficient  improvement  over  the  EKF  with  bias  correction  terms  to  warrant  any 
further  investigation  for  reference  system  applications.  Any  further  examination  of 
second-order  filter  techniques  for  pseudolites  could  be  of  value  for  indoor  pseudolite 
applications,  where  the  nonlinearities  are  significant  enough  to  warrant  such  high- 
order  nonlinear  filtering. 

4  -1  Optimal  Smoothing 

Optimal  smoothers  not  only  increase  the  performance  of  the  floating-point  fil¬ 
ter  and  the  ability  to  resolve  ambiguities,  but  they  also  increase  the  true  window 
over  which  ambiguity  techniques  can  be  applied  in  a  post-processing  application. 
Typically,  the  EKF  took  up  to  7  minutes  (half  of  the  simulation)  to  converge  on 
the  floating-point  solution.  That  is  why  ambiguity  resolution  techniques  were  only 
applied  in  the  second  half  of  the  test  run.  Optimal  smoothing  enables  ambiguity 
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resolution  at  the  start  of  the  mission  if  real-time  processing  is  not  of  critical  impor¬ 
tance. 

The  section  on  optimal  smoothers  in  Chapter  3  detailed  exactly  how  the 
smoothing  algorithm  was  modified  for  base  prn  changes,  psendolites  coming  into 
view,  and  psendolites  going  out  of  view.  In  order  to  test  the  smoother  for  a  double 
difference  base  PRN  change,  the  base  was  manually  changed  from  3  to  8  and  the 
visibility  plot  is  shown  in  Figure  4.18. 
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Figure  4.18  Visibility  Plot  for  20  Kilometer  Maximum  Range 

The  smoother  allows  a  converged  solution  right  at  the  start  of  the  test  mission 
and  can  be  seen  in  Figure  4.19.  This  figure  shows  the  RMS  ambiguity  error  in  the 
top  subplot  and  the  RMS  standard  deviations  of  the  ambiguities  in  the  bottom  plot. 
The  forward  filter  is  depicted  with  a  solid  line,  while  the  smoothed  estimates  are 
shown  with  a  dotted  line. 
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GPS  Week  Seconds-Local  Time 

Figure  4.19  Filter  and  Smoother  RMS  Ambiguity  Error  and  Ambiguity  l-/sigma 
Plots 

These  plots  demonstrate  the  the  smoother  could  handle  the  pseudolites  coming 
into  and  going  out  of  view  along  with  the  changing  base  prn  for  double  differenc¬ 
ing.  The  spikes  in  the  plots  are  a  result  of  a  pseudolite  coming  into  view.  Recall 
that  the  variance  assigned  to  a  pseudolite  that  just  came  into  view  is  rather  large 
in  comparison  to  the  remaining  variances,  which  affects  the  RMS  of  the  standard 
deviations. 

4-8  Summary 

This  chapter  first  described  the  simulation  parameters,  scenario  descriptions, 
and  test  case  definitions.  This  background  was  required  to  set  the  stage  for  the 
single  run  and  Monte  Carlo  results.  In  the  single  run  analyses,  the  filter  was  eval¬ 
uated  based  on  the  position  and  ambiguity  estimation  errors  in  addition  to  using 
the  flight  vehicle  velocities  and  accelerations  to  explain  some  of  the  results  from 
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the  residuals.  The  Monte  Carlo  simulations  evaluated  a  single  frequency  filter  ver¬ 
sus  one  using  a  widelane  observable.  The  Monte  Carlo  simulations  also  included 
a  sensitivity  analysis  of  each  error  source  and  a  comparative  analysis  of  five  filter 
enhancements.  Alternate  flight  trajectories  were  investigated  with  nonlinear  filtering 
techniques.  Lastly,  the  optimal  smoother  was  shown  to  increase  the  epochs  available 
for  ambiguity  resolution. 


4-36 


V.  Conclusions  and  Recommendations 


5. 1  Overview 

This  research  presented  the  theory,  models,  and  simulation  results  for  a  pseudolite- 
based  flight  reference  system.  Previous  research  has  indicated  that  pseudolites  can 
be  used  successfully  for  positioning  and  ambiguity  resolution.  This  research  concen¬ 
trated  on  the  application  and  adaption  of  GPS  carrier-phase  differential  techniques  to 
pseudolite  measurements  for  a  flight  reference  system  application.  The  adaptations 
were  required  due  to  the  differences  in  pseudolite  versus  GPS  navigation. 

The  baseline  algorithm  consisted  of  an  extend  Kalman  filter  that  used  a  dou¬ 
ble  differenced  code  and  carrier-phase  measurement.  Widelane  or  single  frequency 
measurements  could  be  used  in  this  filter,  in  addition  to  a  number  of  possible  fil¬ 
ter  enhancements.  These  filter  enhancements  included  code  averaging  (when  two 
codes  are  available),  bias  correction  terms  (emulating  second  order  filtering),  opti¬ 
mal  smoothing,  and  two  methods  for  reducing  the  residual  tropospheric  delay  that 
exists  after  a  tropospheric  model  has  been  applied.  The  first  method  implemented  a 
weighted  measurement  covariance  matrix  based  on  the  tropospheric-model-predicted 
delays.  The  second  method  for  reducing  residual  tropospheric  error  was  explicitly 
modelling  it  in  the  filter. 

5.2  Conclusions 

A  single  run  of  the  filter  was  evaluated  to  show  typical  performance  of  the 
floating-point  filter.  The  performance  was  investigated  through  analysis  of  the  po¬ 
sition  and  ambiguity  accuracies,  in  addition  to  the  velocities,  accelerations,  and 
residuals.  Although  the  primary  objective  was  to  evaluate  ambiguity  resolution  per¬ 
formance,  position  accuracy  was  also  important.  The  accuracy  of  the  algorithm  in 
the  position  domain  is  also  relevant  if  a  pseudolite-only  system  is  desired  (i.e.,  one 
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based  solely  on  pseudolite  signals).  In  general,  the  floating-point  filter  was  typically 
accurate  to  within  10-15  centimeters  for  the  position  and  two-tenths  of  a  cycle  for  the 
ambiguities  (for  single  frequency  measurements).  For  widelane  measurements,  the 
ambiguity  accuracy  was  within  a  tenth  of  a  cycle,  with  a  position  accuracy  within 
20-25  centimeters.  The  performance  was  investigated  by  single  run  analysis  and 
supported  by  Monte  Carlo  simulations.  For  tested  trajectories,  the  filter  required 
5-7  minutes  to  converge.  Converged  estimates  were  available  at  the  start  of  the  run 
when  an  optimal  smoother  was  applied. 

The  influence  of  the  code  measurements  and  the  relative  motion  between  trans¬ 
mitters  and  receivers  are  the  two  factors  that  allow  the  filter  to  converge  on  the 
ambiguities.  When  the  code  was  essentially  removed,  by  increasing  the  measure¬ 
ment  covariance  values  to  an  extremely  high  value,  the  filter  was  forced  to  rely  on 
the  relative  motion  only.  This  actually  increased  the  speed  of  convergence  with  no 
degradation  of  accuracy.  This  suggested  that  relative  motion  was  the  primary  driver 
for  ambiguities’  convergence.  This  was  confirmed  when  the  filter  processed  measure¬ 
ments  from  a  stationary  trajectory,  which  decreased  the  accuracy  of  the  ambiguity 
estimates  by  an  order  of  magnitude.  This  suggests  that  a  pseudolite-based  flight 
reference  system  does  not  require  code  measurements,  although  code  measurements 
will  add  robustness  to  the  system. 

Monte  Carlo  simulations  were  also  conducted  to  evaluate  filter  performance 
further.  This  analysis  included  an  evaluation  of  widelane  versus  single  frequency, 
a  sensitivity  analysis  of  each  error  source,  and  a  comparative  analysis  of  five  filter 
enhancements.  The  widelane  measurement  reduced  the  magnitudes  of  the  ambi¬ 
guity  errors,  at  the  expense  of  increasing  the  errors  in  the  position  domain.  As  a 
result,  ambiguity  resolution  was  more  easily  conducted  with  a  widelane  frequency 
than  with  a  single  frequency  implementation.  The  filter,  using  widelane  frequency 
measurements,  was  able  to  resolve  100  percent  of  the  ambiguities  correctly,  while  the 
filter  using  single  frequency  measurements  was  able  to  resolve  96.3  percent  correctly 


5-2 


and  .5  percent  incorrectly.  The  difference  was  more  dramatic  for  larger  measure¬ 
ment  error  cases.  Given  correctly  resolved  ambiguities,  a  single  frequency  solution 
has  smaller  measurement  errors  than  a  widelane  frequency  solution.  The  filter  per¬ 
formed  well  enough  to  assert  that  a  single  frequency  may  be  all  that  is  required  in  a 
fielded  system. 

The  expected  level  of  measurement  errors  was  increased  and  decreased  by  a 
factor  in  order  to  characterized  filter  performance  to  different  levels  of  measurement 
error.  First,  the  filter  was  not  altered  and  evaluated  for  each  of  the  different  levels 
of  measurement  errors.  This  represented  a  filter  mistiming  (i.e.,  either  the  filter 
was  over-  or  under-estimating  its  performance),  which  allowed  the  robustness  of  the 
filter  to  be  tested.  The  second  type  of  test  involved  re-tuning  the  filter  in  order  to 
evaluate  its  performance  with  the  correct  level  of  measurement  error.  The  re-tuned 
filter  improved  the  filter  performance  from  99.3  percent  resolved  correctly  to  100 
percent  for  the  best  case,  and  59.2  percent  resolved  correctly  to  63.8  percent  for  the 
worst  case. 

The  sensitivity  analysis  for  each  error  source  suggested  that  the  pseudolite  po¬ 
sition  errors  and  the  residual  tropospheric  error  were  the  dominant  error  sources. 
Great  care  should  be  taken  when  surveying  the  antenna  positions  of  the  pseudolites 
and  reference  receiver.  The  filter  sensitivity  to  the  un-modelled  tropospheric  delay 
error  motivated  two  of  the  filter  enhancements  developed  as  part  of  this  research — the 
weighted  measurement  covariance  matrix  and  tropospheric  model  error  states.  The 
weighted  measurement  covariance  matrix  method  utilized  the  tropospheric  model 
output  in  selecting  measurement  uncertainty  values  based  on  the  predicted  tropo¬ 
spheric  delay.  This  relative  weighting  of  measurement  uncertainty  was  based  on  the 
range  and  elevation  angle  (i.e.,  the  longer  the  range  and/or  lower  the  elevation  angle, 
the  larger  the  predicted  tropospheric  delay  and  thus  uncertainty  magnitude). 

The  second  enhancement  explicitly  estimated  the  tropospheric  model  error  as 
an  additional  state  in  the  filter.  The  error,  when  expressed  as  a  percentage,  was  the 
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modelled  as  a  first  order  Gauss  markov  process.  Initially,  the  error  percentage  of 
the  mobile  receiver  and  reference  receiver  were  modelled  separately.  The  filter  could 
successfully  estimate  the  percentage  of  tropospheric  model  error  in  the  absence  of 
measurement  noise,  multipath,  and  pseudolite  location  errors.  When  these  error 
sources  were  included  in  the  measurement-corrupted  ranges,  the  filter  correctly  esti¬ 
mated  the  mobile  error  percentage,  but  not  the  reference  error  percentage.  The  lack 
of  motion  between  the  reference  receiver  and  the  pseudolites  caused  the  pseudolite 
location  errors  to  be  bias-like.  The  filter  was  apparently  unable  to  distinguish  the 
tropospheric  model  error  from  the  biased  multipath  and  pseudolite  location  error 
present  at  the  reference  receiver.  This  problem  suggested  the  use  of  a  single  tropo¬ 
spheric  state  that  included  the  error  percentage  of  both  receivers.  This  implementa¬ 
tion  successfully  improved  ambiguity  resolution,  and  was  the  only  enhancement  that 
reduced  errors  in  the  position  domain. 

The  other  filter  enhancements  included  code  averaging,  second  order  filter¬ 
ing,  and  optimal  smoothing.  Neither  the  code  averaging  nor  the  bias  correction 
terms  significantly  enhanced  the  ambiguity  resolution  process.  The  bias  correction 
terms  from  a  second  order  filter  were  implemented  to  improve  the  linearization  ap¬ 
proximation  for  the  extended  Kalman  filter.  The  effect  of  the  measurement  model 
nonlinearities  are  attributed  to  two  factors.  The  first  is  the  degree  of  uncertainty 
of  the  receiver,  specifically  in  a  direction  orthogonal  to  the  line  of  sight  between  a 
transmitter  and  receiver.  The  second  factor  associated  with  nonlinearities  is  short 
transmission  ranges  which  result  in  more  spherical  wavefronts  of  the  received  signals. 
Multiple  trajectories  for  a  practical  flight  reference  system  were  tested  to  investigate 
the  impact  of  measurement  model  nonlinearities  and  the  benefit  of  the  addition  of 
bias  correction  terms  or  second  order  filters.  It  was  concluded  the  level  of  uncer¬ 
tainty  in  the  receiver  position  is  so  small,  when  compared  to  the  line-of-sight  ranges, 
that  bias  correction  terms  and  second  order  filters  are  not  helpful  for  the  tested  tra¬ 
jectories.  If  the  uncertainty  is  significant  when  compared  to  the  ranges  (like  would 
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be  the  case  for  indoor  pseudolite  applications),  further  exploration  of  second  order 
nonlinear  filtering  techniques  is  warranted. 

The  optimal  smoothing  algorithm  developed  by  Meditch  [20]  is  not  tolerant 
to  the  state  and  covariance  matrices  changing  dimensions  between  measurement 
updates.  The  number  of  ambiguity  states  varies  between  epochs  based  on  the  number 
of  visible  pseudolites.  This  is  caused  by  pseudolites  coming  into  view  and  going  out  of 
view.  The  existing  smoothing  algorithm  cannot  handle  the  state  dimension  changes. 
The  smoother  calculates  a  state  vector  and  covariance  matrix  at  the  current  time 
ti,  based  on  the  filter  state  and  covariance  at  the  current  time  tll  and  both  the 
filter  and  smoother  outputs  at  the  future  time  ti+1.  As  part  of  this  research,  the 
smoothing  algorithm  was  modified  to  allow  the  number  of  ambiguity  states  to  change 
between  epochs.  The  smoothing  algorithm  could  not  tolerate  a  state  to  represent  one 
quantity  at  one  epoch  and  a  different  quantity  at  the  next.  For  GPS  and  pseudolite 
applications  this  will  occur  when  a  base  double  difference  PRN  goes  out  of  view.  The 
smoothing  algorithm  was  also  modified  to  allow  a  change  in  the  base  double  difference 
PRN.  A  translation  matrix  was  formed  based  on  the  base  PRNs  at  two  adjoining 
epochs  and  the  vector  of  visible  pseudolites.  The  translation  matrix  was  applied  to 
both  the  state  vector  and  covariance  matrix  to  allow  the  smoothing  algorithm  to 
form  the  smoothed  state  and  covariance  estimate  properly. 

The  position  and  ambiguity  solutions  from  the  floating-point  filter  suggest  that 
the  pseudolite  ambiguities  can  be  resolved  with  a  pseudolite-only  system,  and  that 
further  integration  with  an  INS  or  other  measurement  sources  is  not  required  to 
obtain  high-accuracy  position.  This  research  concluded  that  carrier-phase  measure¬ 
ments  with  resolved  ambiguities  can  be  produced  from  pseudolite  signals  to  incor¬ 
porate  and  improve  accuracy,  especially  during  periods  of  GPS  jamming,  to  flight 
reference  systems. 
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5.3  Contributions 


This  thesis  has  provided  contributions  in  various  areas  of  research.  The  fol¬ 
lowing  list  details  these  contributions. 

•  The  largest  contribution  of  this  thesis  is  the  development  and  testing  of  the 
navigation  filter  for  a  pseudolite-based  flight  reference  system.  Although  very 
similar  to  a  GPS  navigation  algorithm,  this  design  concentrated  on  the  differ¬ 
ences  between  GPS  and  pseudolite  systems. 

•  The  Meditch  optimal  smoothing  algorithm  was  modified  to  handle  a  state  and 
covariance  matrix  with  changing  dimensions,  ft  was  also  adapted  for  changes 
in  the  quantities  that  are  represented  in  the  matrices.  These  adaptations  can 
be  easily  extended  to  GPS  navigation  applications. 

•  The  weighted  measurement  covariance  matrix  was  developed  as  part  of  this 
research.  This  method  takes  advantage  of  tropospheric  model  predictions  to 
weight  the  expected  variance  of  the  measurements.  This  significantly  improved 
the  ambiguity  resolution  of  the  filter. 

•  The  tropospheric  model  error  states  were  also  developed  as  an  alternative 
method  for  reducing  tropospheric  delay  error.  This  included  an  additional 
state  that  also  aided  ambiguity  resolution.  This  enhancement  also  improved 
the  accuracy  of  the  filter  in  the  position  domain. 

•  This  research  included  the  analysis  of  the  bias  correction  terms  and  two  second 
order  filters  to  reduce  the  effect  of  measurement  model  nonlinearities.  The 
analysis  indicated  periods  when  divergence  can  occur  with  bias  correction  terms 
and  ways  to  eliminate  the  divergence. 

5-4  Recommendations 

The  filter  in  this  research  performed  well  in  its  ability  to  resolve  pseudolite 
carrier-phase  ambiguities,  and  further  research  and  evaluation  for  this  concept  is 
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warranted.  The  following  recommendations  are  included  to  provided  to  guide  the 
next  logical  steps  for  a  pseudolite-based  flight  reference  system. 

•  Apply  the  filter  to  real  pseudolite  flight  test  data.  It  will  be  necessary  to  find 
a  suitable  reference  system  for  comparison.  Pseudolite  data  is  available  from 
the  original  Holloman  test  along  with  tape  measure  results  [25].  However, 
the  geographic  size  of  this  test  was  fairly  small,  and  it  will  therefore  not  fully 
challenge  the  filter  with  significant  tropospheric  error. 

•  Incorporate  an  inertial  system  to  compare  the  performance  of  a  Pseudolite /Inertial 
Navigation  System  (PL/INS)  to  a  pseudolite-only  system.  The  coupling  of  an 
inertial  system  with  pseudolite  measurements  should  only  improve  the  perfor¬ 
mance  of  the  floating-point  filter,  and  would  be  helpful  if  there  are  pseudolite 
measurement  dropouts.  The  Kalman  filter  used  in  this  research  would  have 

to  be  converted  to  an  error-state  Kalman  filter  that  estimates  the  error  in  the 
inertial  system. 

This  recommendation  is  relatively  easy  to  implement  if  real  flight  test  data  is 
available  for  both  GPS  and  an  INS.  The  fixed-integer  solution  from  the  GPS 
receiver  could  be  used  as  the  true  trajectory  in  order  to  simulate  pseudolite 
measurements  that  are  corrupted  with  simulated  measurement  errors.  The 
filter  would  use  the  real  INS  data  and  the  simulated  pseudolite  measurements. 

•  Implement  a  Multiple  Model  Adaptive  Estimation  (MMAE)  algorithm  for  the 
set  determination  function  of  the  ambiguity  resolution  process.  This  would 
replace  the  ratio  test  of  the  residuals  with  parallel  Kalman  filter  conditioned 
on  each  possible  set  of  ambiguities  [15]. 

•  MMAE  techniques  can  also  be  applied  to  reducing  the  residual  tropospheric 
error.  This  could  be  implemented  instead  of  the  weighted  measurement  covari¬ 
ance  or  the  tropospheric  model  error  states. 
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•  Model  each  visible  pseudolite  with  a  separate  state  for  the  tropospheric  model 
error  state  method.  If  the  time  constant  is  lowered,  the  filter  may  be  able  to 
reduce  the  pseudolite  position  error  and  multipath. 

•  Develop  a  practical  method  that  selects  the  best  period  of  time  to  perform  the 
ambiguity  search  process.  This  decision  should  be  based  on  pseudolite  visibil¬ 
ity,  magnitude  of  acceleration  of  test  vehicle,  and  size  of  ambiguity  covariance 
values.  This  research  applied  ambiguity  resolution  techniques  to  epochs  scat¬ 
tered  throughout  the  entire  data  set,  with  the  exception  of  the  first  few  minutes 
to  allow  the  filter  to  converge.  Selecting  the  best  time  conditioned  on  a  high 
number  of  visible  pseudolites,  small  accelerations  of  the  test  vehicle,  and  low 
ambiguity  covariance  values  should  provide  a  better  methodology  for  a  prac¬ 
tical  system.  If  not  all  the  pseudolites  were  visible  during  this  period,  the 
filter  can  use  the  fixed  ambiguity  solutions  with  reduced  covariances  to  reflect 
that  these  are  the  correct  ambiguities  and  iteratively  solve  for  the  remaining 
ambiguities. 
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